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I   INTRODUCTION 

A.  STATEMENT  OF  PROBLEM 

The  seismic  response  of  shorty  single  or  multiple  span,  highway 
bridges  are  greatly  affected  by  the  phenomenon  of  soil-structure  inter- 
action.  The  dynamic  forces  exerted  by  backfills  on  the  abutments  often 
add  significantly  to  the  maximum  seismic  forces  developed  in  the  overall 
structural  system.  Also,  bridges  of  this  type  usually  have  relatively 
short  and  stiff  columns  which  interact  strongly  with  their  supporting 
foundations.  Neglecting  these  interaction  effects  can  lead  to  large 
errors  in  predicting  design  loads.   It  is  therefore  the  objective  of 
this  investigation  to  establish  appropriate  mathematical  models  which 
will  yield  realistic  seismic  response  for  certain  soil-structure  sys- 
tems of  this  type.  Further,  computer  programs  are  written  to  carry  out 
time-history  dynamic  analyses  as  an  aid  to  the  design  process. 

B.  REVIEW  OF  PERTINENT  RESEARCH  INVESTIGATIONS 

Numerous  analytical  investigations  have  been  carried  out  in  the 
past  to  determine  dynamic  earth  pressures  acting  on  retaining  walls. 
One  of  the  earliest  was  carried  out  by  Okabe  and  later  a  similar  study 
was  reported  by  Mononobe  and  Matsuo  [80,  74]*.   In  each  of  these  inves- 
tigations an  equivalent  static  earth  pressure  was  determined  as  a 
function  of  acceleration  with  its  resultant  acting  on  the  wall  at  a 

♦Numbers  in  square  brackets  refer  to  Bibliography  numbers. 
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position  corresponding  to  that  of  static  fluid  pressure,  i.e.  at  a  height 
one-third  the  distance  from  the  base.  Much  later,  Valera  used  the 
finite  element  method  to  predict  the  intensities  of  seismic  earth  pres- 
sures acting  on  rigid  walls  [106] .  In  this  investigation,  nonlinear 
soil  behavior  was  considered.  Seed  and  Whitman  summarized  31  investi- 
gations carried  out  from  1926  to  1969  and  made  valuable  suggestions 
for  their  application  to  design  [97] .  Later,  Wood  derived  analytical 
expressions  for  earth  pressures  acting  against  rigid  and  rotating  walls 
using  linear  soil  properties  [117] .  All  of  these  investigations  were 
concerned  with  intensities  of  dynamic  pressures  but  in  each  case  they 
assumed  the  pressures  to  act  independently  of  the  dynamic  response  of 
the  wall  itself. 

In  recent  years,  the  influence  of  soil  conditions  on  the  seismic 
response  of  structural  systems  has  received  considerable  attention. 
Seed  recently  published  a  report  covering  two  main  aspects  of  this  prob- 
lem:  (1)  changes  in  seismic  ground  motions  adjacent  to  buildings  as  a 
result  of  physical  interaction  effects,  and  (2)  changes  in  seismic 
response  of  buildings  as  a  result  of  the  changes  in  ground  motions  due 
to  different  soil  deposits.  This  investigation  did  not,  however,  in- 
clude full  dynamic  soil-structure  interaction  effects  for  structures 
other  than  buildings  [98] . 

Few  analytical  investigations  have  been  reported  in  which  dynamic 
soil-bridge  interaction  effects  were  treated  rigorously.  One  such 
study  was  reported  by  Penzien  in  1970  [82] .   In  this  particular  inves- 
tigation, the  soil  foundation  was  represented  by  a  series  of  lumped 
masses  interconnected  by  bilinear  hysteretic  shear  springs  having 
properties  which  varied  with  soil  depth  and  also  interconnected  with 
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viscous  linear  dashpots.   The  bridge  deck,  supporting  piers,  and  pile 
foundations  were  also  modelled  as  lumped  mass  systems.   The  three  dimen-  . 
sional  effects  of  the  foundation  soils  were  determined  using  the  Mind- 
lin  Theory  of  the  elastic  half-space.  Equations  of  motion  were  developed 
which  considered  all  dynamic  interaction  effects.  Recently,  Tseng  and 
Penzien  reported  an  investigation  on  the  seismic  response  of  long"  multi- 
ple span  bridges  in  which  the  mathematical  modelling  included  soil- 
structure  interaction  effects  at  the  bases  of  columns;  however,  due  to 
the  long  structural  types  considered,  soil-abutment  interaction  effects 
were  neglected  [105] . 

As  far  as  experimental  results  are  concerned,  most  investigations 
have  been  carried  out  on  simple  wall-soil  systems.   Important  contribu- 
tions have  been  made  by  Matsuo,  Ishii  and  Arai  and  Tsuchida,  Matsuo 
and  Ohara,  Tschebotanoff ,  Ohara,  and  Niwa  168,  52,  69,  104,  79,  78], 

Shepherd  and  Charleson,  and  Shepherd  and  Sidwell  have  conducted 
field  experiments  on  an  existing  bridge  in  the  small  loading  range 
[100,  101].  Their  results  show  that  the  soil  around  the  bridge  displayed 
high  energy  absorption  qualities,  i.e.  provided  considerable  damping  to 
the  overall  soil-structure  system.  Broms  and  Ingelson  have  also  con- 
ducted field  experiments  to  determine  the  static  earth  pressures  acting 
on  abutment  walls  [11] . 

C.   REVIEW  OF  DAMAGES  CAUSED  BY  EARTHQUAKES 

Iwasaki,  Penzien,  and  Clough  prepared  an  extensive  summary  of  the 
damages  caused  to  bridges  during  earthquakes  for  the  period  1923  to 
1971  [53] .  Jennings  and  Wood,  and  Lew,  Legendecker,  and  Dikkers  have 
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reported  on  the  damages  to  bridges  during  the  San  Fernando  earthquake 
of  February  9,  1971  [59,  65] .   Duke  and  Leeds  reported  similar  damages 
caused  by  the  Chilean  earthquake  of  I960,  while  Rose,  Seed,  and  Migliac- 
cio  reported  on  the  Alaskan  earthquake 'of  1964  [29,  92] .  As  indicated 
in  these  reports,  short  bridges  have  suffered  damages  ranging  from 
column  failures  to  cracking,  tilting,  and  even  overturning  of  piers, 
abutments,  and  wing  walls.  Clearly,  this  evidence  shows  that  large 
dynamic  forces  are  induced  in  the  overall  structural  system  by  backfill 
earth  pressures. 

D.   SCOPE  OF  PRESENT  INVESTIGATION 

The  present  investigation  is  a  study  of  the  seismic  response  of 
short,  single  or  multiple  span,  highway  bridges  having  narrow  abutment 
backfills  as  shown  in  Fig.  1. 

Since  the  earth  pressures  acting  against  the  abutments  can  greatly 
affect  the  seismic  response  of  the  bridge,  the  mathematical  modelling 
includes  a  two-dimensional  soil  element  representation  of  the  abutment 
backfills.   This  representation  accounts  for  nonlinear  soil  properties 
and  allows  different  vertical  boundaries  to  be  present  as  shown  in 
Figs.  ID,  IE  and  IF.  The  bridge  deck,  piers  (or  columns),  and  abut- 
ments are  modelled  using  prismatic  beam  elements  which  may  be  permitted 
to  have  hysteretic  yielding  properties.  A  frictional  element  is  used  to 
model  the  nonlinear,  discontinuous  behavior  at  the  interfaces  of  back- 
fill soils  and  abutments  and  nonlinear,  discontinuous  type  of  expansion 
joint  element  is  included.  The  soil  foundation  flexibilities  under  the 
columns  are  represented  by  equivalent  columns.  The  mathematical  model 
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of  the  overall  soil-structure  system  permits  the  study  of  interaction 
effects  and  yields  the  distribution  of  forces  throughout  the  system. 
Chapter  II  of  this  report  describes  the  different  elements  used 
in  the  mathematical  modelling  and  presents  the  derivations  of  elastic 
stiffnesses  used  for  these  elements.  Chapter  III  discusses  soil  mate- 
rial properties  used  in  the  modelling  with  particular  emphasis  on  non- 
linear properties.  Chapter  IV  develops  the  stiffnesses  of  nonlinear 
elements  with  emphasis  on  the  derivation  of  the  elasto-plastic  force- 
displacement  relations  for  the  soil  finite  elements  and  the  bridge 
column  elements.   Chapter  V  presents  the  numerical  techniques  used  * 
the  time-history  dynamic  analysis  while  Chapter  VI  presents  the  r   it 
of  parameter  studies.  Chapter  VII  presents  certain  conclusions  and  re- 
commendations. Finally,  listings  of  the  computer  program  are  presented 
in  Appendix  A. 
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II   ELASTIC  STIFFNESSES  OF  MODEL  ELEMENTS 

Six  basic  elements  are  used  in  the  mathematical  modelling  of  the 
bridge-soil  system  CD  soil  finite  element  of  a  continuum,  (2)  soil 
boundary  element,  C3)  frictional  element  at  interface  of  soil  and  abut- 
ment wall,  (4)  prismatic  beam  element,  (5)  expansion  joint  element, 
and  (6)  equivalent  column  of  foundation.  It  is  the  purpose  of  this 
chapter  to  describe  these  elements  and  to  derive  their  elastic  stiff- 
nesses. 

A«   SOIL  FINITE  ELEMENT  OF  A  CONTINUUM 

The  bridge-soil  system  shown  in  Fig.  1  is  subjected  to  earthquake 
motions  in  the  vertical  direction  and  in  one  horizontal  direction  coin- 
ciding with  the  longitudinal  axis  of  the  bridge.  Because  of  the  condi- 
tions of  symmetry,  it  is  assumed  that  the  soils  adjacent  to  the  abut- 
ments respond  to  these  motions  in  a  two-dimensional  manner.  This 
assumption  allows  these  soils  to  be  modelled  using  two-dimensional  fi- 
nite elements  of  a  continuum  interconnected  at  their  nodal  points.  Ma- 
terial properties  are  defined  individually  for  each  element  which  may 
have  an  arbitrary  quadrilateral  or  triangular  shape. 

The  quadrilateral  element  used  is  the  isoparametric  element  for 
which  the  geometry  and  displacements  are  described  in  terms  of  the 
same  parameters  of  similar  order.  Using  the  natural  coordinate  system 
and  interpolation  displacement  models,  the  isoparametric  formulation  has 
several  advantages  over  the  generalized  coordinate  method  [25].  First, 
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the  nodal  displacements  provide  a  more  direct  visualization  of  actual 
structural  deformations  than  do  the  generalized  displacements.  Second, 
it  is  computationally  more  efficient  since  no  transformation  from  one 
coordinate  system  to  another  is  required.  Finally,  the  local  coordinate 
system  of  the  isoparametric  approach  coincides  with  the-  global  coordi- 
nate system;  thus,  eliminating  the  need  for  transformation  of  loads  and 
stiffnesses  from  one  coordinate  frame  to  the  other. 

The  following  derivation  of  the  stiffness  matrix  for  the  two- 
dimensional  quadrilateral  isoparametric  element  is  similar  to  that  found 
in  certain  textbooks  125,  34,  118] . 

1.  Coordinate  System  -  The  global  Cartesian  coordinates  (x,y)  and 

the  local  (or  natural)  coordinates  (s,t)  as  shown  in  Fig.  2  are  related 

through  an  interpolation  function  h.  by  the  relations 

k 
x  =  E  h.  x. 

k 

y  =  E  h.  y. 

where 


h  =  H   (1-s)  (1-t) 

l 

h.  -  h   (1+s)  (1-t) 
2  (2) 

h  -  h   (1+s)  (1+t) 

3 

h  =  h   (1-s)  (1+t) 


2.   Strain-Displacement  Equations  -  The  displacements  are  approxi- 
mated using  the  same  interpolation  function  which  can  provide  for  both 
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flexible  and  rigid  body  modes.   In  this  case,  one  can  state 


u   (sft)  =  .£,  h.  u 
x  1=1   i   x 


u   (s,t)  =  .£.,  h.  u 
Y         i=l   i  YL 


(3) 


The  two-dimensional  strains  are  obtained  by  taking  derivations  of  the 
displacements  with  respect  to  x  and  y  in  the  following  manner: 


xx 


e  =  /g 


yy 


xy 


(4) 


3h. 
L   ,   1  i      ' 

1=1   dy   x.    dx   y, 


Since  functions  h.  (i=l,2,3,4)  are  expressed  in  terms  of  s  and  t,  the 
chain  rule  of  derivations  can  be  applied  using  the  equations 


oh. 

1 

8x 

ah. 
i 

as 

a       ah. 

OS       1 

ax     at 

3t 
3x 

oh. 

_  l 

3h. 

l 

0       ah. 

OS      _1_ 

3t 

3y 


3s 


3y   3t 


ay 


(5) 


The  chain  rule  is 


3(  )■  ) 
3s   ( 

3x 
3s 

3s 

\  3(  ) 
J   3x 

3(  )  f 

at 

3x 
3t 

3y_ 
3t 

)  3(  ) 
(  3y 

(6) 
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which  in  its  inverted  form  becomes 


3(   ) 


3x 
3C  ) 


3v. 

3t 

3s 

f  3(  ) 
\  3s 

1 

) 

J 

3x 
~  3t 

3x 
3s 

)  3<  ) 
/  3t 

3y 


where  the  Jacobian  J  is  defined  by 


j  =  a*      h.   _  ?JL  h. 
3s   3t     3t  3s 


Thus,  the  required  derivatives  are  given  by 


3s 
3x 

3t 
3x 

1^ 

3y_ 

3t 

-3v_ 
3s 

3s 
3y 

3t 
3y 

J 

-3x 
3t 

3x 
3s 

Taking  derivatives  of  Eqs.  (1) ,  one  can  write 

3h. 


3x 
3s  " 

^    i 
£  — *—  x. 

ds   l 
1=1 

3x 
3t  " 

f  3h. 

£  -57-  x- 
3t  l 
1=1 

3s 

«♦   3h. 

3s  i 

1=1 


3       «t  3h. 
^  =   I 


3t 


1=1 


3t  *i 


(7) 


(8) 


(9) 


(10) 


Making  use  of  Eqs.  (2)  and  (10),  the  Jacobian,  as  defined  by  Eq.  (8) 
becomes 
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J   =   8  tS(X3  -  V(yi  -  Y2]     '     (X1  ~  V(y3  7!V. 

+  t(x2  -  x3)  (y2  -  y.)  -  (xx  -  xj  (y2  -  y^    (11) 

+   (xx  -  x3)  (y2  -  yj  -  (x2  -  Xlf)  (yx  -  y3)]  - 

and  the  strain  vector  defined  by  Eq.  (4)  becomes 


£  = 


3n. 
l 

3x 

0 

3h 

2 

3x 

0 

3h 

3 

3x 

0 

3h 
<3x 

0 

0 

3h 
i 

3y 

0 

3h 

2 

3y 

0 

3h 

3 

3y 

0 

3h 
3y 

3h 
l 

3h 
l 

3h 

2 

3h 

2 

3h 

3 

3h 

3 

3h 

3h 

3y 

3x 

17 

3x 

3y 

3x 

3y 

3x 

or  in  compact  form  may  be  expressed  as 


IX 


u 


lYi 


u 


2X 


U 


2Y 


u 


3X 


U 


3y 

i 
^x 

1 


(12) 


e  =  B  u 


(13) 


The  derivatives  in  the  coefficient  matrix  of  Eq.  (12)  can  now  be 
expressed  in  their  final  maniputable'  form 


3h 
l 

3x 
3h 

2 

3x 


3h 


3x 

3h 
«»_ 

3x 

3h 
l 

3y 


3       _ 


~-  I  y   (h  +  h)  +y    (h-h)  -  y   (h  +  h)] 

16J  212  3        «f  2  «f       1  h 

r^-  [-y    (h  +  h  )  +  y   (h  +  h  )  +  y    (h -  -  h  )  ] 

lbJ  112  323  <fl3 


r^-r  [y    (h-h)    -y    (h  +  h)   +y    (h  +  h)] 

16J  12  h  22  3«»3  «f 


y—  [y    (h  +  h)  +y    (h-h)  -  y    (h  +  h  )] 

16J  11  i»  23  1  3       «♦  3 

t^t  [-x    (h  +  h  )  +  x    (h  -  h  )  +  x    (h  +  h  )  ] 

16J  212  3        2  «f  «♦        1  «t 
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8h 

-57-     =     TTT  I   x    (h  +  h)    -x    (h  +  h)    +  x    (h-h)] 

dy  16J  112  323  4        3  1 

-T1    -     TZT  I  x    Ch  -  h  j   +  x   (h  +  h  )    -  x    (h  +  h  )] 

dy  16J  1       if  2  2       2  3  4       3  4 

3h 

-5^-    =     — -  [-x    (h  +  h  )    +  x    (h  -  h  )    +  x    (h  +  h  )] 

dy  16J  114  213  334 

(14) 

Using  standard  index  notation  with  i  permutating  1  through  4,  Eqs-  (14) 
can  be  written  in  the  two  single  equations 


dh. 

-£■  =  tIt  I  y-^  <h-  +  h-A  )  +  y-j.  Ch..  -  h.A  )  -  y.A  (h.+  h.A  1] 

dX       16J    'l+l   1    1+1       1+2   1+3    1+1       1+3   1    1+3 

3h.      , 

-T  =  TFT  I-y-j.  <h-  +  h-u.  >  +  *-A  0».A  -  h_  )  +  x_  (h.+  h.^  )] 

dy      16J    ■'l+l   1    1+1       1+2   1+1    1+3       1+3   1    1+3 

(15) 


3.  Element  Stiffness  and  Numerical  Integration  -  Considering  a 
thickness  03  (normal  to  the  plane  of  element  shown  in  Fig.  2)  ,  the  ele- 
ment stiffness  matrix  can  be  expressed  in  the  form 

K  =  /  B_TC_  B  dv        /   B_TC_  B_  OJdA  (16) 

vol  area 

where  C_  is  the  stress-strain  matrix  with  integration  being  carried  out 
over  the  entire  area  of  the  element.  For  purposes  of  numerical  integra- 
tion, Eq.  (16)  can  be  written  in  terms  of  the  s  and  t  coordinates  giving 

1  1     T 
K  =   /  /  B  C  B  (det  J)  Uids   dt  (17)  . 

-l  ^l 
Upon  application  of  standard  one-dimensional  numerical  integration 
formulas,  Eq.  (17)  becomes 
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K  =  E  E  a),  a),   (det  J)  B  (slf  t)  CB  (s.f  t)w    (18) 
-     -k   ]  k         -   3   k 3   k 

in  which  s .  and  t  are  integration"  points  and  o) .  and  to  are  the  appro- 
priate weighting  functions.  Using  the  Gauss-Legendre  quadrature  formula 
of  degree  2,  one  obtains 


s.  =  +  0.577350 
l   — 

t.  =  +  0.577350 

l   — 

03.  =  1 

1 


i  =  1,2 


(19) 


For  plane  stress,  matrix  C  has  the  form 


£=r^ 


1   V    0 

V 

0 


1   0 

l-v 


(20) 


while  for  plane  strain  it  becomes 


-  "   (1+V) (1-2V) 


1-V   V 
V   1-V 
0    0 


0 
1-2V 


(21) 


where  E  is  Young's  Modulus  and  V  is  Poisson's  ratio. 

As  shown  above,  the  final  form  of  the  element  stiffness  matrix  is 
a  function  of  geometry,  element  thickness,  and  material  properties. 
The  overall  stiffness  matrix  of  the  entire  system  is  assembled  by  the 
direct  stiffness  method  [22] .   The  band  width  of  this  matrix  depends 
upon  the  system  of  numbering  nodal  points.   Therefore,  the  nodal  points 
should  be  numbered  in  a  manner  to  minimize  computer  storage  requirement 
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and  operating  time. 

B.   SOIL  BOUNDARY  ELEMENT 

Boundary  elements  are  used  in  the  mathematical  model  to  account 
for  the  elastic  action  which  occurs  at  the  vertical  boundaries  of  the 
soils  being  considered.   The  type  of  element  used  for  this  purpose  is  a 
simple  elastic  spring  as  shown  in  Fig.  3„  This  element  has  a  stiffness 
matrix  k  in  the  local  coordinate  system  of  the  form 


k  =  E 


-1 


-1 


(22) 


where  E  is  the  equivalent  stiffness  of  the  spring,  which  can  be  obtained 
using  standard  methods  [82,89] . 

C.   PRISMATIC  BEAM  ELEMENT 

Prismatic  beam  elements  are  used  to  model  the  bridge  deck,  piers, 
abutment  walls  and  equivalent  columns  for  soil  foundations.  The  deri- 
vation of  the  stiffness  matrix  for  a  beam  element  considering  axial , 
shear,  and  bending  deformations  can  be  found  in  many  textbooks  [86] . 
Therefore,  only  the  main  features  of  this  matrix  will  be  presented  here. 

The  force  components  acting  on  the  beam  element  are  axial  forces 
s  and  s  ,  shearing  forces  s  and  s  ,  and  bending  moments  s  and  s  as 

1         H  2         5  3d 

shown  in  Fig.  4.   The  positive  directions  of  these  force  components  (s) 
and  their  respective  displacement  components  (u1)  correspond  to  those 
shown  in  the  figure.   The  stiffness  matrix  for  this  element  in  terms  of 
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the   local  coordinates   is  given  by  the   standard  form 


—(6X6) 


EA 
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12EI 


A3  (l+cf    ) 


6EI 


z2a+i  ) 
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EA 
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-12EI 
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Ty        z 
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£2(1+<|>    ) 
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(2-d    )EI 
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A (1+4    ) 
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SYMMETRICAL 


EA 


12EI 


£3 (1++  ) 

y 


-6EI 


£2(l+<!>  ) 

y 


(4+<j>  )EI 

y     ; 

£(14  ) 
y 


where 


12EI 


G  Asy  £' 


24U+V)  ~-  (-4) 
Asy    £ 


C23) 


(24) 


represents  the  shear  deformation  parameter  for  reinforced  concrete 
elements  and  where 


£  = 


z 

A  = 

Asy  = 


element  length 

radius  of  gyration  about  z  axis 

Moment  of  inertia  about  z  axis 

cross  sectional  area 

equivalent  shear  area  in  y  direction 
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E  =  modulus  of  elasticity 

G  =  modulus  of  elasticity  in  shear 

V  =  Poisson's  ratio 

The  matrix  equation  relating  displacements  in  the  local  coordinates 
u*  to  displacements  in  the  global  coordinates  u_  is  given  by 


u' 


X  u 


(25) 


where  X  has  the  two-dimensional  form 


X  = 


^l 

y* 

0 

0 

0 

0 

V* 

Dx/D* 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

V** 

y      I 

0 

0 

0 

0 

-y0* 

y    * 

0 

0 

0 

0 

0 

0 

l 

(26) 


and  where 


D  = 


D  = 


D„  = 


x. 

D 

-  x. 

1 

yj 

-yi 

^D 

2    +   D 

2 

(27) 


Thus,  the  stiffness  matrix  K  in  global  coordinates  becomes 
K  =  XT  k  X 


(28) 


as  indicated  above,  the  stiffness  matrix  is  linearly  proportional  to 
the  areas  and  moments  of  inertia.   Therefore,  a  unit  width  analysis  can 
be  carried  out,  if  desired,  by  dividing  all  areas  and  moments  of 
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inertia  by  the  width  of  the  bridge  deck. 

D.   FRICTIONAL  ELEMENT 

A  so  called  frictional  element  is  used  to  model  the  frictional 

action,  separation,  and  impact  which  take  place  at  the  interfaces  of 

* 

soil  backfills  and  abutment  walls.  ,  This  element  has  the  following 

characteristics  (1)  the  frictional  force  per  unit  area  is  proportional 

to  the  normal  interface  pressure  and  a  coefficient  of  friction;  thus, 

slippage  occurs  when  the  direction  angle  of  the  resultant  of  pressure 

and  friction  exceeds  the  soil  angle  of  friction,  (2)  impact  occurs  at 

the  interface  upon  closure  of  any  gap  which  may  have  earlier  developed, 

and  (3)  no  frictional  resistance  can  develop  at  the  interface  when  wall 

and  soil  surfaces  have  separated.  Discontinuous  elements  similar  to 

this  have  been  developed  by  Ghaboussi  and  Wilson,  Scholes  and  Strover, 

White  and  Enderly,  and  Tseng  and  Penzien  [40,  94,  111,  105];  however, 

the  element  developed  by  Goodman  and  Taylor  has  been  adopted  here  [41] . 

This  frictional  element  is  a  four  nodal  element  as  shown  in  Fig.  5 
having  length  L  but  with  a  height  equal  to  zero,  i.e.  nodal  points  1 
and  4  coincide  as  do  points  2  and  3.   It  is  shown  in  a  local  coordinate 
system  with  the  origin  at  the  center  of  the  element  and  the  x  axis  di- 
rected along  the  length  of  the  element. 

The  relative  displacement  vector  u_  is  expressed  in  terms  of  the 
displacement  vector  u.   through  the  linear  interpolation  function  for- 
mula 
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u  = 
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(29) 


The  material  property  matrix  k_  expressing  the  stiffness  per 
unit  length  in  the  normal  and  tangential  directions  is  given  by 


k  = 


k     0 

s 


n 


(30) 


Upon  applying  the  variational  principle  of  solid  mechanics,  the  elastic 
stiffness  matrix  in  the  local  coordinate  system  becomes 
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(31) 
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Ec   EXPANSION  JOINT  ELEMENT 

An  expansion  joint  may  be  present  between  the  bridge  deck  and  an 
abutment  as  shown  in  Fig.  6.  This  joint  can  develop  horizontal  fric- 
tional  forces  which  should  be  modelled  properly.  For  this  purpose, 
the  frictional  element  previously  described  can  be  adopted;  however, 
two  additional  boundary  conditions  must  be  imposed,  namely,  (1)  the 
relative  displacement  of  the  upper  and  lower  part  of  the  joint  element, 
u   -  u  ,  cannot  cause  a  gap  closure  (between  deck  and  abutment)  great- 

3X     2X  - 

er  than  the  original  gap  dimension,  and  (2)  frictional  forces  can  be 
developed  only  when  the  relative  displacement,  u   -  u  ,  causing  a 

IX     3  X 

widening  of  the  gap  is  less  than  the  original  "seat"  dimension.  As  soon 
as  the  relative  displacement  (u   -  u  )  exceeds  the  original  seat  di- 
mension, the  bridge  deck  falls  from  its  support;  thus,  the  computer 
analysis  is  stopped  at  this  point. 


F.   EQUIVALENT  COLUMN  FOR  FOUNDATION 

Various  mathematical  models  have  been  used  for  structural  founda- 
tions. As  reported  in  the  literature,  Parmelee,  Whitman  and  Roesset, 
Dobry,  and  others  [81,  108]  used  a  simple  spring  dashpot  model;  Jennings 
and  Bielak  and  Richart,  Hall,  and  Woods  [58,  89]  used  the  elastic  half 
space;  Whitman  [109]  used  an  equivalent  lumped  mass  model;  and,  Dans 
and  Butterfield,  Finn,  Lysmer  and  Kuhlemeyer,  and  Wilson  {31,  36,  66, 
114]  used  a  finite  element  mesh.  Various  methods  have  been  reported 
for  estimating  the  lateral  stiffness  of  foundations  employing  piles 
[84,  88] . 
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The  replacement  of  the  foundation  flexibility  with  an  equivalent 
column  through  use  of  the  elastic  half  space  theory  and  under  the 
assumption  of  quasi-static  behavior  has  been  reported  by  Penzien,  et. 
al.  [8] .  This  method  has  been  adopted  in  the  present  investigations. 

The  first  step  in  finding  the  equivalent  column  for  the  founda- 
tion is  to  determine  the  lateral,  vertical,  and  rotational  stiffnesses 
of  the  foundation  at  the  footing  (or  pile  cap)  level.  Once  these  three 
stiffnesses  have  been  determined,  the  foundation  is  replaced  by  a  column 
of  length  L,  flexural  stiffness  EI,  and  axial  stiffness  AE  which  when 
fixed  at  its  base  provides  the  equivalent  lateral,  vertical,  and  ro- 
tational stiffnesses  to  the  footing.  The  stiffness  matrix  for  this 
column  is  given  in  the  form  of  Eq.  (23)  under  the  assumption  of  no 
shear  deformation,  i.e.  <J>  =0. 

y 

The  foundation  stiffnesses  can  be  obtained  by  either  the  numerical 
procedure  outlined  by  Penzien  or  by  a  closed  form  approach  reported  by 
Gerrand  and  Harrison  [82,  39] . 

1.  Closed  Form  Solution  fpr  Lateral  Stiffness  of  Circular  Foot- 
ing -  The  closed  form  solution  for  lateral  stiffness  is  available  for 
a  circular  shaped  footing  as  shown  in  Fig.  7-c.  If  the  footing  is 
rectangular  in  shape,  an  equivalent  radius  r  should  be  calculated  for 
a  circular  footing  having  the  same  area,  then 
8  r  G 

\  ■  sk-  (32> 

where 

k   =  lateral  stiffness  for  a  single  footing 
r   »  radius  of  footing 
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G  =  shear  modulus  of  soil 
V  =  Poisson's  ratio 

For  a  pier  of  multiple  supports,  the  interactions  between  footings  are 
estimated  by  Eq.  (33)  which  describes  the  lateral  displacement  at  a 
distance  r  from  a  loaded' footing. 


1   r,„  ..,   .-11.  „    (r2-l) 


u  =  T>,  TZ  \    t(2-V)-sin   -  +  V  -  ■   >    '    )  C33) 

h  4TrrQG  r         r2 

where 

u  =  the  displacement 

T,  =  the  total  applied  force  at  a  single  footing 
n 

r  =  the  distance  in  terms  of  radius  r 

o 

By  calculating  the  average  displacement  of  all  footings  in  a  pier  and 

dividing  the  total  force  by  the  average  displacement,  the  equivalent 

stiffness  of  a  pier  with  multiple,  footings  is  then  obtained. 


2.   Numerical  Solution  Using  the  Mindlin  Equation  to  Calculate 
Lateral  Stiffness  -  The  Mindlin  equation  has  the  following  form. 
When  the  x  component  of  displacement  on  the  surface  is  to  be  calculated, 
as  produced  by  a  single  concentrated  force  P  located  at  the  origin 
(o,  o,  o)  on  the  surface  of  an  isotropic  half  space  and  acting  in  the 
x-direction 


u  (x,y,o)  .  Pte,o.o)    U-M  +L+4g-v)(i-2v) 

x      '*'  16lTCl-V)G       ]        R  R  R 

+  x2    [   I     +      3~4v   _  4 (1-V) (1-2V) 
R3             R3  R3 
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(34) 


in  which 

R2   =   x2  +  y2 

Under  the  assumption  that  the  lateral  force  is  uniformly  distributed 
over  the  area  of  the  footing,  the  procedure  to  calculate  the  lateral 
stiffness  of  a  single  footing  is  (1)  replace  the  uniform  pressure  by 
as  many  concentrated  loads  P  as  accuracy  requires,  (2)  calculate  the 
lateral  displacement  of  each  load  point  over  the  footing  area  due  to 
each  concentrated  load,  (3)  sum  the  resulting  displacements  at  each 
load  point  within  the  footing  as  caused  by  the  full  set  of  loads,  (4) 
average  the  load  point  displacements,  and  (5)  divide  the  total  resultant 
force  by  the  average  of  the  load  point  displacements  to  obtain  the 
lateral  stiffness. 

Again,  if  a  pier  has  many  footings,  the  same  averaging  procedure 
is  applied  by  assuming  a  concentrated  force  at  the  center  of  each  foot- 
ing; thus,  the  lateral  stiffness  of  each  pier  foundation  can  be  obtained. 

3.   Comparison  of  Two  Methods  for  Lateral  Stiffness  -  To  compare 
the  two  previously  described  methods  for  obtaining  foundation  lateral 
stiffness,  consider  the  pier  shown  in  Pig.  8  having  four  footings. 
Assuming  the  footing  is  8.5  x  8.5  ft  and  the  distance  between  footings 
is  26  ft,  the  individual  footing  stiffness  k  and  pier  stiffness  k 
obtained  by  the  two  approaches  have  the  following  values  in  k/ft: 

!32.5  G     Numerical  method 
25.5  G     closed  form  solution 

f 100.6  G    Numerical  method 

P   1 

/  80.0  G     closed  form  solution 
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2 

Where  the  unit  of  G  is  in  k/ft  . 

In  each  case  the  stiffnesses  differ  by  approximately  20%.   In  the 

numerical  solution  the  footing  is  considered  to  have  36  concentrated 

forces  applied  at  equal  grid  intervals  over  the  area.   If  the  resultant 

force  had  been  discritized  at  more  than  36  points,  the  differences  in 

» 
stiffness  would  have  decreased. 

4.  Vertical  Stiffness  with  Friction  Piles  -  Assuming  the  friction 
force  per  unit  length  of  pile  varies  in  a  linear  manner  from  a  maximum 
value  at  the  top  to  zero  at  the  bottom,  and  assuming  zero  vertical  dis- 
placement at  the  bottom  of  the  pile,  the  vertical  stiffness  of  the  pile 


k  is 
v 


k   =  3AE/L  (35) 


Where  L,  A,  and  E  are  the  length,  area,  and  modulus  of  elasticity, 
respectively,  of  the  pile. 

5.  Vertical  Stiffness  Without  Piles  -  Assuming  uniform  vertical 
displacements  as  shown  in  Fig.  7-a,  the  vertical  stiffness  of  a  circular 
footing  without  piles  is  given  by 

4  r  G 

k   =   .  °  (36) 

v     1-v 

For  other  footing  shapes,  similar  stiffness  relations  have  been  reported 
by  Lysrrer  and  Duncan  [67]  . 

6.  Rotational  Stiffness  Without  Piles  -  Assuming  rotational 
displacements  as  shown  in  Fig.  7-b,  the  rotational  stiffness  of  a  cir- 
cular footing  without  piles  is  given  by 
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8  r  3  G 


3(1-V) 


(37) 


The  rotational  stiffness  of  a  rigid  footing  with  piles  can  be  calculated 
directly  from  the  vertical  stiffness  of  each  pile. 
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(a)  UNIFORM  VERTICAL  DISPLACEMENT 
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Ill   MATERIAL  PROPERTIES  OF  SOIL 

In  this  chapter,  the  pertinent  soil  material  properties  used  to 
establish  a  non-linear  finite  element  model  are  discussed.   The  non- 
linear stress-strain  relations,  Mohr's  envelope  and  p-q  diagram,  the 
concepts  of  active  and  passive  stresses,  and  the  damping  characteristics 
of  soil  are  described. 

A.   STRESS-STRAIN  RELATIONSHIP 

Basic  to  establishing  the  force-displacement  relationships  for 
soil  elements  as  employed  in  the  mathematical  model  are  the  stress- 
strain  relationships  for  the  materials  involved.  These  materials  in- 
clude both  cohesionless  and  cohesive  soils. 

Numerous  authors,  including  Bishop  and  Henkel,  Bishop,  Comforth 
and  Seed,  have  reported  procedures  for  determining  the  shear  modulus  of 
sand  [9,  10,  23,  96] .   The  investigators  show  that  the  shear  modulus 
is  strongly  influenced  by  confining  pressure,  strain  amplitude,  and 
void  ratio.   In  the  present  investigation,  the  equivalent  secant  shear 
modulus  as  determined  by  extreme  points  on  the  hysteresis  loop,  Fig.  9, 
is  adopted.   This  modulus  can  be  estimated  using  the  relation  proposed 
by  Seed  [96] ,  namely, 

.  \ 
G  =  1000  k   (O*  )  psf  (38) 

2   m  c 

where  G  is  the  secant  modulus  of  sand,  0*   is  the  effective  mean  stress 

m 

and  k  is  a  parameter  which  depends  upon  void  ratio  e  and  strain  amplitude 
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as  shown  in  Fig.  10.   Triaxial  tests  show  that  k  depends  only  upon 

2 

void  ratio  e  at  very  low  strain  levels  (y  ^  10  '  percent) .  At  inter- 
mediate strain  levels  (10   <^  y   <_  10   percent)  ,  it  still  depends  pri- 
marily upon  void  ratio  but  is  also  slightly  influenced  by  state  of 
stress  and  the  friction  angle  <}>  .   At  very  high  strain  levels,  k  is 

2 

essentially  a  constant  which  is  almost  independent  of  state  of  stress, 
friction  angle,  and  the  void  ratio.  Thus  for  practical  purposes,  the 
values  of  k  can  be  assumed  to  vary  only  with  strain  amplitude  and  void 
ratio  as  shown  in  Fig.  10. 

For  prescribed  values  of  void  ratio 'and  confining  pressure,  the 
stress-strain  relationship  is  non-linear  with  the  shear  modulus  chang- 
ing with  shear  strain.  Upon  the  initiation  of  yielding  the  modulus 
becomes  very  small.  A  tangent  modulus  curve  transformed  from  the  equi- 
valent shear  modulus  curve,  Fig.  10,  is  shown  in  Fig.  12.   Theoretically, 
one  could  use  a  revised  tangent  modulus  over  each  time  interval  of  a 
dynamic  analysis;  however,  it  is  believed  to  be  more  practical  to  use 
a  more  simplified  form  of  stress-strain  relationship.   In  the  present 
investigation  a  trilinear  stress-strain  relationship,  Fig.  13-a,has  been 
adopted  with  the  shear  modulus  remaining  constant  during  each  of  three  dif- 
ferent loading  stages.   In  the  initial  stage,  the  soil  element  is  in 
its  geostatic  state,  i.e.  the  vertical  stress  is  equal  to  the  weight 
of- soil  (per  unit  area)  above  the  point  of  interest  as  given  by 

y 

<J      =  /  w  dy  (39) 

o 

where  y  is  the  depth  and  w  is  the  unit  weight  of  soil.   The  lateral 
stress  in  the  principal  horizontal  direction  is 
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a   =  k  a  (40.) 

h      o   v 

while  in  the  orthogonal  horizontal  direction  (z)  the  stress  O     for 
plane  stress  is 


0       =  0  (41) 

z 

and  for  plane  strain  is 

k 

a      =  u(a  +  a. )   =  .  °  (a  +ka)  =  k  a  =  ex    (42) 
z       v    h    l+k    v    o  v       o  v    h 

o 

where  ]i  is  Poisson's  ratio.   The  shear  modulus  in  this  stage  is  evalua- 
ted by  Eq.  (38)  with  k  selected  in  accordance  with  the  very  low  strain 
level  shown  in  Fig.  10,   In  the  second  stage  the  soils  in  the  backfill 
are  no  longer  in  the  geostatic  state  due  to  loadings  from  the  bridge. 
In  this  case  the  shear  modulus  is  calculated  using  Eq.  (38) ,  but  k 
is  revised  to  be  consistent  with  the  maximum  shear  strain  in  the  ele- 
ment.  The  third  stage  occurs  after  the  initiation  of  yielding  in  which 
case  the  stress-strain  relation  is  evaluated  in  accordance  with  a  theory 
of  plasticity  as  described  in  Chapter  IV. 

In  those  cases  where  curves  of  k  vs.  Y  as  shown  in  Fig.  10  are 

2 

unavailable,  a  bilinear  stress-strain  relationship  is  adopted  with 
the  initial  modulus  being  estimated  by  an  empirical  method  and  the 
second  stage  modulus  being  calculated  by  the  theory  of  plasticity  as 
shown  in  Fig.  13-b. 

Turning  our  attention  now  to  certain  clay  materials,  test  data 
show  that  at  very  low  strain  levels,  the  shear  modulus  varies  almost 
linearly  with  shear  strength  [112] .   In  a  summary  report,  Seed  has 
presented  a  curve  of  secant  shear  modulus  versus  shear  strain  for 

-31- 


saturated  clays  {96] .   The  shear  modulus  expressed  by  this  curve,  Fig. 

14 ,  is  in  the  normalized  form  G/s  where  s   is  the  undrained  shear 

u        u 

strength  of  the  material.  A  curve  showing  the  degradation  of  secant 
shear  modulus  with  shear  strain,  as  obtained  by  a  number  of  investiga- 
tors, is  shown  in  Fig.  11  [96].  These  relationships  can  be  used  as  a 
guide  for  estimating  the  initial  modulus  of  clay  when  laboratory  test 
data  on  the  specific  material  are  unavailable. 

Fig.  15  gives  a  tangent  modulus  curve  for  clay  which  is  consistent 
with  the  secant  modulus  curve  of  Fig.  11.  Test  results  show  that  the 
shear  modulus  is  essentially  independent  of  the  confining  pressure. 
A  re-evaluation  of  shear  modulus  before  the  initiation  of  yielding  be- 
comes much  less  important  for  this  material  than  for  sand.  For  this 
reason,  a  bilinear  stress-strain  relation  has  been  assumed  for  cohesive 
soils  in  the  present  investigation.  This  initial  modulus  is  estimated 
from  the  data  previously  described  and  the  tangent  stiffness  after 
yielding  is  obtained  using  the  same  basic  procedure  as  used  for  sands. 
No  hardening  effects  after  the  initiation  of  yielding  are  considered. 

B.   MOHR'S  ENVELOPE,  p-q  DIAGRAM 

The  most  widely  used  yield  criterion  in  soil  mechanics  is  the 
Mohr -Coulomb  criterion  which  relates  the  normal  stress  O     and  the 
shear  stress  T  at  the  failure  plane  by  the  equation 

X   =  c  +  0     tan  <)>  (43) 

where  c  is  the  cohesion  of  the  soil  and  ty   is  its  friction  angle.   The 
soil  coefficients  are  usually  obtained  by  conducting  triaxial  tests  at 
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various  confining  pressures  and  by  drawing  a  common  tangent  line  to 

the  resulting  Mohr's  circles  as  shown  in  Fig.  16.   it  should  be  noted 

that  since  the  confining  pressure  -of  the  triaxial  test  is  uniform, 

i.e.  0     =  0    ,  the  Mohr's  circle  can  be  used  in  its  familiar  two- 
2     3 

dimensional  form.   Equation  (43)  when  plotted  on  the  O-T  plane  is  called 
Mohr's  yield  envelope  whidh  implies  (1)  elastic  behavior  for  a  state  of 
stress  whose  Mohr's  circle  lies  entirely  below  the  envelope,  (2)  yield- 
ing, or  impending  yielding,  for  a  state  of  stress  whose  Mohr's  circle 
has  the  envelope  as  its  tangent,  and  (3)  that  any  state  of  stress  whose 
Mohr's  circle  crosses  the  envelope "is  not  permitted. 

An  alternate  way  of  plotting  the  results  of  triaxial  tests  is  to 
plot  the  stresses  at  the  plane  of  maximum  shear  on  the  p-q  plane, 
namely, 

a    +  a 

p  =    1  =   3  (44) 


and 


a     -  a 

q  .=    -  n      3  (45) 


In  Fig.  17,  points  A,  B,  and  C  in  the  0-T  plane  represent  the  stresses 
on  planes  of  maximum  shear.   A  line  drawn  through  these  same  points 
in  the  p-q  plane,  as  shown  in  Fig.  18,  is  called  the  k  line.   The 
equation  of  this  line  is 

q   =  a  +  p  tan  a  (46) 

Strength  parameters  c  and  <f)  can  be  computed  from  Eq.  (46)  through 
the  relations 

sin  (f>  =  tan  a  (47) 
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c  =  a/cos  <f>  (48) 

Co   ACTIVE  AND  PASSIVE  STRESS 

Using  the  previously  defined  yield  criterion,  active  and  passive 
soil  stresses  can  be  defined  as  they  relate  to  the  backfill  pressures 
exerted  on  abutment  walls.   Starting  with  soil  equilibrium  in  its 
geostatic  state  which  has  principal  axes  in  the  vertical  and  horizontal 
directions,  decrease  the  horizontal  compressive  stress  continually  until 
the  shear  strength  of  the  soil  is  reached  and  failure  occurs.  The  hori- 
zontal compressive  stress  at  the  point  of  failure  is  active  stress.  On 
the  other  hand,  if  the  horizontal  compressive  stress  is  continually 
increased,  the  shear  strength  of  the  soil  will  again  be  reached  at  a 
much  higher  stress  level  called  the  passive  stress.   The  Mohr's  circles 
representing  these  two  failure  conditions  are  shown  in  Fig.  19  where 
G  and  a  represent  the  vertical  and  horizontal  stresses,  respectively, 
in  the  geostatic  state  arid  Q  and  0     represent  the  active  and  passive 
stresses,  respectively.   It  should  be  noted  here  that  since  the  two 
horizontal  stresses  are  assumed  equal,  only  two  stress  components  are 
required  to  describe  the  stress  conditions. 

The  stress  conditions  in  the  C-T  diagram  of  Fig.  19  are  again 
shown  in  the  p-q  diagram  of  Fig.  20,  where  point  A  is  the  geostatic 
state  of  stress,  point  B  is  the  passive  state  of  stress,  and  point  C 
is  the  active  state  of  stress.   Lines  AC  and  AB  are  stress  paths  which 
depict  the  successive  states  of  stress  which  occur  in  changing  from 
the  geostatic  condition  to  the  active  and  passive  conditions,  respective- 
ly. 
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When  applying  the  concepts  of  active  and  passive  stresses  to  more 
complex  stress  conditions,  certain  modifications  are  needed.  For  ex- 
ample, the  stress  conditions  of  the  bridge  backfill-soil "treated  in  this 
investigation  is  essentially  three-dimensional.   In  this  case  shear 
stresses  may  exist  on  both  vertical  and  horizontal  planes  even  under 
gravity  loadings.  Referring  to  Figs.  1-A  and  B  for  the  element  coordi- 
nates, the  Mohr's  circle  for  this  condition  may  be  as  shown  in  Fig.  21. 
Further,  it  should  be  noted  that  if  the  side  slopes  of  the  embankment 
as  seen  in  Fig.  1-B  are  not  sufficiently  flat,  the  yield  plane  will 
show  as  a  line  in  the  y-z  plane.   In  the  present  study  however  it  is 
assumed  that  these  slopes  are  sufficiently  flat  so  that  the  yield  plane 
shows  as  a  line  in  the  x-y  plane.   Thus  the  state  of  stress  at  failure 
as  shown  on  an  element  of  soil  in  the  x-y  plane,  can  be  represented  by 
a  Mohr's  circle  tangent  to  the  Mohr's  envelope  as  shown  by  Circle  1 
in  Fig.  22.   States  of  stress  in  the  x-z  and  y-z  planes  would  be  repre- 
sented by  Mohr's  circles  falling  entirely  within  the  above  defined  cir- 
cle as  shown  by  Circle  II  and  III  in  Fig.  22.   Because  of  these  restric- 
tions, the  problem  can  be  treated  as  a  two-dimensional  problem  in  the 
x-y  plane.   Finally,  in  the  present  investigation,  all  stresses  vary 
prior  to  reaching  a  state  of  yield.   Based  on  the  previous  assumptions, 
active  and  passive  states  can  be  defined  in  an  equivalent  manner. 
Referring  to  Fig.  23  where  the  absolute  values  of  the  p  and  q  stresses 
are  plotted,  p  values  at  static  loading  and  failure  conditions  can  be 
compared.   If  |p[  is  assumed  to  decrease  continually  with  increasing 
|q|  stress  until  the  Mohr ' s  envelope  is  reached,  active  failure  occurs. 
On  the  other  hand,  if  [p(  is  assumed  to  increase  continually  with  in- 
creasing |q|  stress  until  the  envelope  is  reached,  passive  failure 
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occurs.   Stresses  p  and  p  represent  the  active  and  passive  stresses 
in  this  case. 


Do   DAMPING 

Soils,  like  all  materials,  exhibit  energy  dissipation  when  subjected 
to  cyclic  loading.  Different  methods  have  been  used  for  measuring  damp- 
ing depending  upon  the  strain  levels  involved  [96]t.  Forced  vibration 

_2 

tests  have  been  used  for  strain  levels  in  the  range  10   to  10   percent, 
free  vibration  tests  have  been  used  for  Strain  levels  in  the  range  10  3 
to  1.0  percent,  and  static  tests  have  been  used  to  measure  the  hystere— 
tic  energy  absorption  for  strain  levels  in  the  range  10  ';  to  5  percent. 
The  energy  dissipation  in  the  latter  case  has  been  expressed  in  terms 
of  equivalent  viscous  damping  ratios  [102,  103] . 

Experimental  evidence  shows  that  two  major  factors  influence  the 
amount  of  damping  exhibited  by  sands,  namely,  confining  pressure  and 
strain  amplitude.  Damping  in  this  case  tends  to  decrease  with  an  in- 
crease in  confining  pressure  and  tends  to  increase  with  strain  ampli- 
tude. Clay  materials,  on  the  other  hand,  exhibit  damping  which  is 
essentially  independent  of  confining  pressure  but,  like  sands,  the 
damping  increases  with  strain  amplitude. 

In  modelling  the  damping  as  measured  in  soils,  it  has  generally 
been  carried  out  by  using  an  equivalent  viscous  damping  system.  An 
important  basic  study  carried  out  along  these  lines,  using  a  wide 
range  of  yield  values,  has  been  reported  by  Hudson  147];  see  Fig.  24. 
In  this  case,  the  equivalent  viscous  damping  ratio  corresponding  to 
elasto-plastic  hysteretic  damping  for  a  single  degree  of  freedom  system 
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is  given. 

In  the  present  investigation,  hysteretic  damping  has  not  been 
converted  to  an  equivalent  viscous'  form.  Rather,  it  is  treated  in  its 
true  strain  dependent  form.   Since  hysteretic  damping  treated  in  this 
manner  accounts  for  energy  dissipation  only  for  strain  levels  above 
yield,  it  is  necessary  to  also  include  viscous  damping  to  represent  the 
low  strain  level  velocity  dependent  damping.  For  this  purpose,  Rayleigh 
damping  has  been  used  in  the  present  investigation. 
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^Strain 


FIG.  9        HYSTERETIC  STRESS -STRAIN  RELATIONSHIPS  AT 
DIFFERENT  STRAIN  AMPLITUDES -SECANT  MODULUS 
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Shear  Strain -percent 
FIG.  10        SECANT  SHEAR  MODULI  OF  SANDS  AT  DIFFERENT  VOID  RATIOS 

(After  Seed  and  Idriss) 
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FIG.  13      IDEALIZED   STRESS -STRAIN  RELATIONSHIP 
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A 


MOHR'S  ENVELOPE 


FIG.  16   MOHR'S  ENVELOP  FOR  SOIL  SAMPLE 
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COMPRESSION 


FIG.  17  STRESSES  AT  MAXIMUM  SHEAR  PLANE 
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FIG.  19    MOHR  CIRCLE 
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Kf-LINE 


Kf-LINE 


FIG.  20    STRESS  PATHS  FOR  RANKINE  ACTIVE 
AND  PASSIVE  CONDITIONS 
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FIG.  21      STRESS  CONDITIONS  AT  FIELD 
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FIG.  22   MOHR  CIRCLE  IN  3  PLANES 


Ps  =  STATIC  STATUS 
Pa  =  ACTIVE  FAILURE 
Pp  =  PASSIVE  FAILURE 
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FIG.  23     ACTIVE   AND  PASSIVE  FAILURE 
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FIG.24     EQUIVALENT  VISCOUS  FRICTION  VS  YIELD  RATIO 
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IV   NON-LINEAR  STIFFNESSES  OF  MODEL  ELEMENTS 

During  periods  of  low  amplitude  oscillation,  a  bridge-soil  system 
can  be  modelled  using  the  linear  elements  as  described  in  Chapter  II. 
It  may  be  necessary  in  this  case  to  treat  the  friction  and  expansion 
joint  elements  in  a  piece-wise  linear  fashion.  For  a  severe  earthquake 
however  inelastic  deformations  may  occur  in  the  concrete  columns  and/or 
backfill  soils,  separations  and  impacts  may  develop  between  the  abutments 
and  backfills,  slippages  may  take  place  in  the  expansion  joints,  and 
yielding  may  take  place  at  the  soil  boundaries  or  in  the  foundation  to 
complicate  the  behavior. 

It  is  the  purpose  of  this  chapter  to  describe  the  non-linear  be- 
haviors of  all  elements  and  to  derive  their  non-linear  stiffnesses. 

A.   GENERAL  ELASTIC -PERFECTLY  PLASTIC  STRESS-STRAIN  RELATIONS 

Presently,  extensive  literature . exists  on  the  theory  of  plasticity 
and  its  application  to  different  types  of  materials  and  structural  ele- 
ments [30,  46,  85],   Recently,  its  application  has  been  extended  to  soil 
structures  and  frame  structures  as  reported  by  Dibaj  and  Penzien,  and 
by  Porter  and  Powell  128,  83]. 

The  first  step  in  deriving  the  stress-strain  relations  for  an  elas- 
tic-perfectly  plastic  material  is  to  assume  a  yield  function  expressed 
in  terms  of  the  stress  space.  This  stress  function  can  be  expressed  as 


f  CT.  .)  =  0  (49) 

ID 
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By  application  of  the  flow  rule,  the  plastic  strain  increment  tensor  is 
derivable  from  this  function  using  the  relation 

6  eij  =  x  8T77  (50) 

3-D 
where  X  is  a  non-negative  scale  factor.  The  total  strain  increment  ten- 
sor is  thus  decomposed  into  its  elastic  and  plastic  components  as  ex- 
pressed by 

Se(.  =  6  e?.  +  6  e*.  (51) 

ID        ID      ID 

The  generalized  Hooke's  law,  relating  the  increment  of  stress  tensor  to 
the  increment  of  elastic  strain  tensor,  can  be  written  in  the  form 

St..  =  C^„  .      6  eE.  (52) 

Using  Eq.  (51) ,  Eq.  (52)  can  be  rewritten  as 

For  an  elastic-perfectly  plastic  material 

6  f  =  || —  6  T .  .   =  0  (54} 

oT . .      IT 

Substituting  Eqs.  (53)  and  (50)  into  Eq.  (54) ,  one  obtains 
3f   -E 


3t.  .  Cijk£ 
ID 


«  S* ■- «  St'  ■  °  l55) 


or 


^z — C...0   6  e  .  -  w- — C...0   kt X  =  0        (56) 

dT.  .      lnkJo  k%  dT .  .      i-jkx      dT,  „ 

lj         J  ID  k^ 
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Solving  for  A  gives 


x  -  hcijkji  !fr:6£k*  <57> 

ID 


where 


1_  E  3f 9f 

h  ■  c«u  3t..   Tr£  <58> 

Substituting  Eq.  (57)  into  Eq.  (50)  results  in  the  relation 

6£ij     =     hCmnkJl     8F-6ekil     8777  (59) 

mn         xj 

A  further  substitution  of  this  relation  into  Eq.  (53)  gives 

St..     =    cE  <6ev0-hCE.         I| — 1| —  6  e     )      (60) 

13  i]kJ6  kJl  xjmn     9t..    9t   -  mn 

or 

6  x..     =    cE.,0      (6  ev0  -  a.       6  e     )  (61) 

lj  xjkx,  kx,         k&mn         mn 

where 

A,  „  =     hC.„  „      7; r (62) 

kJlmn  ijkJl     3t  .  .   3t 

J      x]   mn 

The  final  form  of  the  stress-strain  relation  for  elas to -perfectly  plas- 
tic material  can  now  be  written  in  the  form 


4  Ti3  "  cijk* 6  £kji  (63> 


where 


C  =     C  —     C  A  (64) 

ijk£  ijk£  ijmn     kilmn 
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and  where 


P  E 

C..l0      =     C..        A,  (65) 

ljkA  ljmn     kicmn 


B.   SOIL  FINITE  ELEMENT 

* 

•The  Mohr-Coulomb  criterion  as  stated  in  Eq.  (43)  of  Chapter  m 
may  be  expressed  in  terms  of  the  principal  stresses,  namely 


(a     -  0  )    +    (0     +  a  )    sin  <{>     =     2c  cos  A  (66) 

13  13 

where  0     >   a  >  0     with  tension  as  positive.   In  the  case  of  two 
1—2—3  * 

dimensional  stress  in  the  x,y  plane,  this  relationship  takes  the  form 


{O     +  a   )    sin  <{>  +  2R  -  2c  cos  A  =  0  (67) 

x    y 


where 


R  =  I  (  X  "  Y)2  +  •  (Txy)2]^  (68> 


is  the  radius  of  the  failure  stress  circle  as  shown  in  Fig.  25. 

1.   Tangent  Stiffness  In  Plain  Strain  -  The  derivative  of  tangent 
stiffness  given  by  Eqs.  (49)  through  (65)  will  now  be  presented  in 
matrix  form  for  the  case  of  plain  strain.  Using  the  Mohr-Coulomb  cri- 
terion as  the  yield  function  f(T. .),  and  the  index  notation  for  stresses 

ID 

such  that  T   =  a  ,  T   =  a  ,  and  T   =  T   ,  one  obtains 
11    x   22    y       12    xy 

T   -  t    2         i. 

f  =  (T   +  T   )  sin  <f)  +  2  [  (-^ — — )   +  T2  ]'2  -  2c  cos  6  =  0 

11      22  2  12 

(69) 
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and 


q 


3-D 


3f_ 

8t 


(70) 


ID 


Equation  (50)  becomes 


{SeP} 


X{q} 


(71) 


where 


{q>" 


<   q  >      =     <  q        q        q        > 
^  ^11        22       12 


T      -  T 
<    [  —  sin  <p  +    ( - )  ]  ,  [  —  sin  <p 


x     -  t 


[2-11]    > 


(72) 


Equation  (51)  in  matrix  form  can  be  written  as 


•  {6   e}  =  {8  eE  +  6   eP} 


(73) 


where 


{6  e}  =  <6  z>   =  <  6  e   6  e   6  e   > 

11      22      12 

{6  eP}T  =  <  6  ep  >  =  <  6  ep  6  eP  6  ep  > 


11      22 


12 


{6  eE}T  =  <  6  £E  >  =  <  6  eK  6  eE  6  eE  > 

11      22     12 

{6  x}  =  [CE]    {<5  eE} 


and  where 


{6t}=<6t>  =  <6t       St       6t       > 

11  22  12 


(74) 


[cE] 


E 


(1+V) (1-2V) 


i-v 

V 

0 

V 

l-V 

0 

0 

0 

1-2V 

2 

(75) 
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A  =  h{q}T  [cE]  {6  e} 


(76) 


where 


1 

h 

=       £,*} 

[CE]     {q} 

E 
1+V 

sin     <J) 
2(l-2v) 

1 
h  = 

E   B 
1+V 

<• 

2        2    R2 


m   sin2  »    1_  +  3_  V 
2(l-2v)    2    2   D2 

R 


(77) 

# 

(78) 
(79) 


Equations  (62)  through  (65)  in  matrix  form  become 


[A]   =  h  {q}  {q}T  [CE] 


{6  T>  =   [C]  {6  e} 


(80) 
(81) 


[C]   =   [CE]  -  [CE]  [A] 

=   [CE]  -  h[CE]  {q}  {qT}  [CE] 

=  tcE]  -  [cp] 


(82) 


Letting 


{Q> 


{q>T  [CE] 


<  0  0  0  > 

1+V  *1  *2  *3 

E  sin  $ 

1+V  *  2(1-2V) 


T   "  T     -a 
11    12   sin  9 


T  -  T    T 

11     22     12 


4R    '  2(1-2V) 


4R 


R 


(83) 


,  the  matrix  [C  ]  can  be  expressed  in  the  form 
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[CP]       =      h[CE]     {q>    {q}T    [CE] 


=      h{Q}    {Q}' 


B(l+V) 


Q  Q 

2  1 

Q  Q 

3  1 


Q  Q 


C  Q 

3     2 


Q  Q 

13 
2     3 


(84) 


Finally,  the  tangent  stiffness  in  explicit  matrix  notation  becomes 


IC]  = 


(1+V)  (1-2V) 


1-V   V 
V 
0     0 


1-V    0 
1-V 


b(1+v) 


!A   V,] 


Q  Q   Q  ' 

2  1     2 


2  3 


3  1     3  2     3 


(85) 


2.   Postulate  for  Application  of  Tangent  Stiffness  to  the  Case  of 
Plane  Stress  -  In  order  to  use  the  previously  derived  tangent  stiffness 
for  the  case  of  plane  stress ,   one  must  make  two  assumptions  as  follows: 
(1)  no  yielding  occurs  in  the  third  direction  and  (2)  the  results  of 
triaxial  tests  are  directly  applicable  to  the  case  of  plain  stress. 

Under  these  assumptions,  the  above  derivations  for  plane  strain  also 

E 
apply  to  the  case  of  plane  stress  except  that  the  matrix  [C  ]  is  changed 

to  the  form 


[CE]   = 


1-V' 


1 

V 
0 


V  0 
1  0 
0   1-V 


(86) 


This  results  in  the  following  changes 
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T       -    T  2  T  2 

{q}T    [CE]    {q}    =  -22j  [  (1+V)    sin2c}>   +    (1-V)     (    1  -  4R   22)       +    (1-V)     {-£)    ] 


1        2EB 


h       1-V2 


(87) 
(88) 


,-..    T       -   T  2  .    T  2 

B  =    [  (1+V)    sin2*  +    (1-V)     (  "4R  22)      +    (1-V)     (-^-)    ]    (89) 


{Q>T   =   {q>T    tCE] 


— —    <  Q    Q    Q    > 
l-v2  2     2     3 

1  T       -   T 

<  ~  (1+V)    sin  <f>  +    (1+V) 


1-V 


2         2 


4R 


T      -  T  T 

j   (1+V)    Sin  <J>    -    (1+V)       "4R   ",     (1-V)    -12.  > 


(90) 


[cp] 


=  h[CE]    {q}    {qT}    [CE] 


=  h{Q}    {Q} 


T 


2B(l-v2) 


Q  Q      Q  Q 

12  13 


QQ       Q  ' 

2     1  2 


2     3 


QQ,     QQ      Q  ' 

3     1  3     2  3 


(91) 


Finally,  one  obtains 


[C]    = 


[cE] 


tcp] 


1-V 


1        V 

V 

0        0 


1        0 
1-V 


2B (1+V) 


12  3 


2     1  2 


2     3 


QQ^     QQ       Q  ' 

3     1  3     2  3 


(92) 
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C.   COLUMN  ELEMENT 

1.   Trilinear  Yield  Surface  of  the  Moment-Axial  Force  Interaction 
Diagram  -  Moment-axial  force  interaction  curves  for  the  most  commonly 
used  sections  of  reinforced  concrete  columns  are  available  in  handbooks 
published  by  the  American  Concrete  Institute  [2] .  These  sections  include 
the  spirally  reinforced,  circular  and  square  columns,  and  the  symmetri- 
cally reinforced,  rectangular  tied  columns.  Three  typical  sections  are 
shown  in  Fig.  27.  For  otner  types  of  sections,  a  computer  program  has 
been  developed  to  obtain  the  interaction  curves  by  a  direct  analysis 
method  [105] . 

There  are  three  controlling  points  on  the  interaction  curve  which 
can  be  used  to  approximate  its  form  using  a  trilinear  relationship. 
These  points  are  the  minimum  eccentricity  point  B,  the  balanced  point  C, 
and  the  pure  moment  point  D,  as  shown  in  Fig.  26.   Segment  AB  having  a 
horizontal  slope  defines  the  ultimate  axial  load  capacity  as  that  axial 
load  given  by  the  ACI  code  for  the  minimum  eccentricity  condition  [35] . 
Segment  BC  defining  the  compression  failure  zone  connects  point  B  with 
the  balanced  point  C  which  corresponding  to  a  concrete  strain  of  0.003 
and  a  steel  strain  equal  to  the  yield  strain.  Finally,  segment  CD 
connects  point  C  with  point  D  which  corresponds  to  the  yield  moment  in 
the  presence  of  no  axial  load.  This  latter  segment  defines  the  tension 
failure  zone.   Since  tension  seldom  occurs  in  the  columns  of  short 
bridges,  it  is  not  necessary  to  define  the  interaction  diagram  in  the 
negative  (tension)  region  of  P.  Line  segment  OD  therefore  is  considered 
the  boundary  line  for  the  yield  surface  which  signifies  zero  tension 
capacity. 
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Using  the  approximate  trilinear  form  of  the  interaction  curve,  the 
normalized  yield  stress  function  can  be  written  as 


f.    (S    ,   S   )    =  a.    S     +  b.    S     +  c.    =  0 
112  11  12  i 


(93) 


where 


s    -2- 

1      p2 


s    =^L 

2'     M2 


(94) 


For  line  segments  AB,  BC,  and  CD,  coefficients  a.,  b.,  and  c.  are, 
respectively , 


M 
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1 
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2 
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"  M      ~  1 
2 
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~  P 
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2 

P,          M 

1  1 

=  P0     "  M 

2  2 

C 
3 

P          M 
_    _3 3 

P      "  M 

2              2 

(95) 


2.   Tangent  Stiffness  -  Using  the  trilinear  interaction  relation- 
ship, the  derivation  of  the  elastic-perfectly  plastic  tangent  stiff- 
ness of  a  column  follows  the  same  procedure  used  previously  for  soils 
except  one  must  consider  that  (1)  plastic  deformations  are  concentrated 
at  the  ends  of  the  element  with  the  deformations  taking  place  indepen- 
dently over  zero  lengths  at  each  end  of  the  element,  (2)  plastic  defor- 
mations are  independent  of  the  shear  forces  present,  and  (3)  the  stiff- 
nesses are  expressed  in  terms  of  element  end  forces  and  displacements 
rather  than  in  terms  of  stress  and  strain  as  in  the  case  of  soil  ele- 
ments.  In  this  case,  one  obtains 
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(duP>  = 


du! 


du. 


{q}1   {0} 


{0}    {q} 


=  tq]  U> 


(96) 


P       P 
where  du  and  du  are  the  plastic  deformation  increments  at  ends  I  and 
I       J 

J  of  the  element,  respectively,  and  X  and  X  are  the  associated  pro- 
portionality  factors.   It  follows  therefore  that 


tdu  j     =  <  du   du   du   du   du   du  > 
l     2     3     if     5     6 


W 


,3f: , 


m 


=<a.  Ob.  >  ^  m  =  I  or  J 
l    l  m 


(97) 
(98) 


m 


{du}  = 


{dS>  = 


du. 


/du. 


dS. 


dS. 


,du~ 


'du. 


[KE] 


(99) 


(100) 


where  [K  ]  is  the  elastic  stiffness  matrix  appearing  in  Eq.  23.   Further, 
one  obtains 


{X}  = 


=      [h]     [q]T    [KE]    {du} 


(101) 


[h]'1      =      [q]T    [KE]     [q] 


[A]  =       [q]     [h]     [q]T    [KE] 


(102) 
(103) 


{dS}        =      [K]     [du] 


(104) 
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[K]   =   [KE]  -  [KE]  [A] 

=   [KE]  -  [KE]  [q]  [h]  [q]T  [KE] 

=   [KE]  -  [KP]  (105) 


3 .  An  Approximation  Used  in  Numerical  Iteration  -  Due  to  the 
occurence  of  impact  upon  the  f rictional  element  at  the  interface  of  soil 
and  abutment  wall  elements,  the  time  interval  used  in  a  dynamic  analysis 
to  obtain  a  stable  numerical  solution  must  be  quite  small.   It  must  be 
sufficiently  small  so  that  "overshooting"  of  the  interaction  diagram 
during  a  single  interval  as  shown  in  Fig.  28  is  minimized.  Even  though 
this  interval  is  kept  small,  the  error  introduced  into  a  solution  by 
"overshooting*  has  been  corrected  [16,  60,  105] . 

In  the  present  investigation,  a  simple  procedure ' has  been  adopted 
for  the  overshooting.   If,  as  shown  in  Fig.  28,  the  elastic  stress 
state  assumed  during  an  interval  moves  the  applied  forces  from  point  A 
to  point  B  then  a  transition  from  an  elastic  to  a  yield  state  is  indi- 
cated. While  the  new  force  vector  S   .   as  represented  by  point  B  is 
adopted  without  correction,  point  C  which  is  the  intersection  of  the 

new  stress  vector  S   .  and  the  yield  segment  EF  is  used  to  calculate 

P 
the  slope  of  the  plastic  deformation  vector  du  .  This  same  procedure 

is  used  when  overshooting  occurs  at  a  discontinuity  point  on  the  inter- 
action curve  as  shown  in  Fig.  29.   In  this  case,  the  elastic  stress 
state  assumed  during  an  interval  moves  the  applied  forces  from  point  A' 
to  point  B'.   The  new  force  vector  S   .   intersects  line  EF;  thus,  the 
slope  of  EF  is  used  to  calculate  the  plastic  deformation  increment. 
Other  investigators  have  used  somewhat  different  procedures  for  this 
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correction  [73,  77,  83,  107]. 


D.  FRICTIONAL  ELEMENT 


The  non-linear  behavior  of  the  frictional  element  is  described  in 

terms  of  normal  and  shear  stiffnesses  k  and  k  during  three  distinct 

n      s      3 

stages,  namely,  (1)  when  reparation  occurs,  in  which  case  k  =  k  =0, 

n    s 

(2)  when  compression  occurs  at  the  interface  but  the  shear  strength  of 
the  element  is  not  exceeded  in  which  case  k  and  k  are  assigned  high 
values,  and  (3)  when  compression  occurs  but  the  shear  strength  of  the 

element  is  exceeded  in  which  case  k  =  0  and  k  retains  a  high  value. 

s         n  ' 

The  shear  yield  strength  of  the  element  can  be  defined  by  the  Mohr- 
Coulomb  yield  criterion,  i.e., 


T  =  a  tan  <J>  (106) 

T  w 

where  <J>  is  the  angle  of  wall  friction  between  the  soil  and  the  abut- 
ment wall. 

Before  discussing  the  value  of  <f>  to  be  used,  two  terms  must  be 

Vr 

defined  [63] .  Firstly,  the  constant  volume  frictional  angle  is  defined 
a  -  O 

—  1     1    3 

by  d)   =  sin    ( )   ,  where  O     and  O     are  the  axial  stress  and 

cv  a    +  O        cv         1       3 

1    3 

confining  pressure,  respectively,  in  a  triaxial  test  at  that  stage  when 
the  sand  strains  without  further  volume  change.   Secondly,  peak  friction 
angle  <J>  is  defined  as  the  slope  of  Mohr  envelope  which  is  a  function  of 
the  stresses  indicated  at  the  peak  of  stress-strain  curve  of  a  triaxial 
test. 

For  backfill  against  a  concrete  wall,  Lamb  suggests  that  the  angle 
of  wall  friction  <b  is  about  equal  to  <J>   and  that  it  typically  has  a 
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numerical  value  of  about  30°  [63] .   Seed  and  Whitman  in  a  discussion  of 
dynamically  active  pressure  against  walls  suggest  that  ({)  =  <|)/2  is 
satisfactory  for  most  practical  purposes  [97] . 


E.   EXPANSION  JOINT  ELEMENT 


The  stiffnesses  of  the  expansion  joint  element  can  be  defined  in 

terms  of  normal  and  shear  stiffnesses  k  and  k  during  three  stages, 

(1)  when  the  frictional  resistance  between  the  deck  and  abutment  is  not 

exceeded  in  which  case  k  and  k  are  assigned  high  values,  (2)  when  the 

n      s 

frictional  resistance  is  exceeded,  but  the  gap  between  deck  and  abut- 
ment as  shown  in  Fig.  6a  is  not  closed,  in  which  case  k  =  0  and  k 

s         n 

retains  a  high  value,  (3)  when  the  frictional  resistance  is  exceeded 
and  the  gap  is  closed  in  which  case,  if  the  relative  displacement  indi- 
cated is  consistent  with  gap  closure,  high  values  k  and  k  are  retained, 


F.   SOIL  BOUNDARY  ELEMENT  AND  EQUIVALENT  COLUMN  OF  FOUNDATION 

The  non-linear  behavior  of  the  soil  boundary  element  and  the  equi- 
valent column  foundation  element  can  be  approximated  by  defining  yield 
levels  using  standard  methods  and  adopting  elasto-plastic  hysteretic 
models.   In  the  present  investigation,  only  elastic  behavior  has  been 
considered. 
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FIG.  26  COLUMN  INTERACTION  CURVE 
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(a)    CIRCULAR   SECTION  WITH  BARS  CIRCULARLY  ARRANGED 
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(b)   SQUARE  SECTIONS  WITH  BAR  CIRCULARLY  ARRANGED 
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tc)   RECTANGULAR   SECTION 
FIG.  27     INTERACTION  DIAGRAM  OF  CONCRETE  COLUMN 
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FIG.  28     SITUATION  AT  OVERSHOOTING 
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FIG.  29     OVERSHOOTING  AT  DISCONTINUITY 
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V   DYNAMIC  ANALYSIS  PROCEDURES 

In  the  subsequent  sections,  the  nonlinear  coupled  equations  of  mo- 
tion for  the  discrete  parameter  soil-structure  system  are  formulated 
and  the  method  used  to  establish  their  associated  stiffness,  mass,  and 
damping  matrices  are  described.  Also  presented  are  the  step-by-step 
integration  techniques  employed  in  their  solution. 

Ac   EQUATIONS  OF  MOTION 

Although  some  previous  investigations  have  considered  spacial  vari- 
ations in  the  earthquake  ground  motions  [27,  32,  43] ,  the  present  inves- 
tigation assumes  identical  motions  at  all  base  points  of  the  soil-struc- 
tural system.  This  assumption  is  considered  reasonable  due  to  the  rela- 
tively short  lengths  of  bridges  being  considered. 

The  coupled  equations  of  motion  at  time  t  for  an  N  degree  of  free- 
dom system  subjected  to  rigid  base  excitation  can  be  expressed  in  the 
matrix  form 

[M]  {u>  +  [CI   {u>  ■+    [Kl   {u>  =  {r>  (107) 

t      t         t    t      t 

where  [M]  is  the  constant  mass  matrix,  [C]  ,  and  [K]   are  the  time 
dependent  damping  and  stiffness  matrices,  respectively,  and  where  {u}  , 
{u}  ,  and  iu}  are  the  nodal  point  displacement,  velocity,  and  accelera- 
tion vectors,  respectively.   The  excitation  force  vector  "Cr).  due  to 
rigid  base  motions  is  given  by  the  relation 

{r>  =  -  [M]  {i}  (  u*   +   uY_  )  (108) 

t  gt       gt 
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where  {i}  is  the  unit  vector  and  u  ^  and  u  ^  are  the  horizontal  and  ver- 

gt      gt 

tical  components  of  ground  acceleration. 


Bo   STIFFNESS  MATRIX 

The  complete  stiffness  matrix  [K]   is  assembled  from  the  individual 
element  stiffness  matrices  using  the  direct  stiffness  method  [22] .  The 
individual  element  stiffnesses  during  the  elastic  and  inelastic  ranges 
in  each  time  interval  are  obtained  by  the  procedures  described  in  Chap- 
ters II  and  IV.  The  complete  stiffness  matrix  takes  on  a  symmetric 
banded  form;  thus,  only  the  diagonal  and  the  off -diagonal  terms  on  one 
side  need  be  stored  in  the  computer. 

Co   MASS  MATRIX 

The  mathematical  model  used  assumes  all  mass  as  concentrated  at 
the  nodal  points.   The  diagonal  mass  matrix  which  results  represents  a 
significant  saving  in  computer  storage  and  computational  time  when  com- 
pared with  similar  requirements  for  the  consistent  mass  matrix  121] . 
One-third  of  the  mass  of  each  triangular  element,  one-fourth  of  the  mass 
of  each  quadrilateral  element,  and  one-half  of  the  mass  of  each  beam 
element  are  lumped  at  their  respective  nodal  points.  No  rotational 
moments  of  inertia  are  assigned  to  these  masses.  The  resulting  mass 
matrix  for  the  complete  soil-structural  system  takes  the  form 


IM]  =  diag  <  M  M  M  >  (109) 

12      n 

where  M.  is  the  mass  associated  with  the  ith  degree  of  freedom  and  n 
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is  the  total  number  of  degrees  of  freedom  present  in  the  system.   The 
static  condensation  procedure  is  used  to  eliminate  the  degree  of  free- 
dom of  zero  rotational  masses  in  the  solution  of  eigenvalues. 

D.   DAMPING  MATRIX 

Various  methods  have  been  used  by  investigators  to  determine  the 
viscous  damping  matrix  corresponding  to  matrix  IC]   in  Eq.  (107)  144]. 
Wilson  and  Penzien  have  described  two  methods  for  evaluting  orthogonal 
damping  matrices  [115] .  The  first  method  relates  the  modal  damping 
ratios  to  the  coefficients  in  the  Caughey  series  form  [15] .  The  second 
method  is  a  direct  approach  which  expresses  the  damping  matrix  as  a  sum 
of  a  series  of  matrices  each  of  which  produces  damping  in  only  one  par- 
ticular mode.  The  second  approach  has  the  advantage  that  prescribed 
damping  ratios  in  all  modes  are  easily  controlled. 

The  Rayleigh  damping  matrix  which  constructs  a  damping  matrix  from 
a  scaled  linear  combination  of  the  mass  and  stiffness  matrices  is  used 
in  the  present  investigation.   This  type  of  damping  matrix  has  the  ad- 
vantage that  it  can  be  calculated  directly  using  the  relation 

[CJt  -  cc[M]  +  3[K]t  (110) 

where  a  and  8  are  scalar  quantities  to  be  prescribed.   By  properly 
selecting  these  scalar  values,  the  damping  ratios  can  be  controlled  in 
two  normal  modes.   It  can  be  shown  that  these  quantities  are  related  to 
the  damping  ratios  (£)    and  circular  frequencies  (u))  of  the  ith  and  jth 
normal  modes  through  the  equations 
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2wio   (£  .to.  -  £  to  ) 
ffl  =    i  j  ^j  i   ^ij'  (111) 

(a).2  -  a).2) 
i     D 

2  (£.0).  -  £.w.) 

3  =  i-J: 3-J- 

(oo.2  -  to.2) 


(112) 


Further,  it  can  be  shown  that  if  a  and  3  satisfy  Eqs.  (Ill)  and  (112), 
the  damping  ratio  in  nth  normal  has  the  value  given  by 

<*  +  3  to  2        .  ..... 

c-  _  n_  (113) 

^n      2  (0 
n 

In  the  present  investigation,  the  numerical  values  of  a  and  3  are  deter- 
mined  by  using  the  initial  elastic  soil-structural  system  and  by  pre- 
scribing the  damping  ratios  of  any  two  modes  of  the  system.  These  quan- 
tities are  then  held  constant  at  these  values  throughout  the  time 
history  of  response  including  those  periods  of  time  when  the  system 
responds  inelastically. 

As  shown  by  Eq.  (107)  ,  the  stiffness  matrix  varies  with  time  due 
to  nonlinear  effects;  therefore,  the  damping  matrix  also  varies  with 
time.  Because  the  stiffnesses  in  the  system  decrease  considerably 
during  periods  of  element  yielding,  the  viscous  damping  present  during 
these  periods  also  decreases.   It  should  be  kept  in  mind,  however,  that 
the  major  sources  of  energy  dissipation  during  these  periods  are  the 
hysteresis  loops  in  the  force-deformation  relations  as  described  in 
Chapter  III. 

E.   STEP-BY -STEP  INTEGRATION  TECHNIQUES 
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Having  the  solution  of  the  coupled  equations  of  motion  at  time 
t,  the  step-by-step  integration  procedure  allows  one  to  obtain  their 
solution  at  a  later  time  t+At.  To'  develop  this  procedure,  the  matrix 
equation  of  motion  is  transformed  to  its  incremental  form  by  subtract- 
ing Eq.  (107)  for  time  t  from  the  corresponding  equation  for  time  t+At 
as  given  by 

[M]     ({u}.    +  {Au}.)    +    [CI       ({u>     +  {Au>   )    +   [Kl       ({u>     +  {Au}j 
t  t  tt  t  tt  t 

=   {R>t  +    {AR>t  (114) 

In  this  equation,  the  incremental  quantities  represent  those  changes 
taking  place  during  the  interval  At  following  time  t.  Thus,  one  ob- 
tains the  incremental  form 

[M]    {Auh    +    [CI      {Au>     +    [Kl      {Au>     =  {Ar>  (115) 

t  t  t  t  t  t 

To  find  the  incremental  changes,  various  procedures  can  be  employed 
[7,  76] .   The  differences  in  these  procedures  relate  to  the  analytical 
form  of  the  variation  in  response  over  the  time  interval  At.   In  the 
present  investigation,  two  different  analytical  forms  have  been  pro- 
grammed for  computer  solution,  namely,  constant  acceleration  and  linear 
acceleration.   These  forms  lead  to  the  following  equations  for  velocity 
and  displacement  at  time  t+At  expressed  in  terms  of  the  state  vectors 
at  time  t  and  the  acceleration  vector  at  time  t+At: 

Constant  Acceleration 

{i}t+At  "  {i}t  +  J  4t  mt  +  \  At  K}t+At  (116) 
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{u)t+at  -  {u}t  +  &t  {i}t  +  \  At2{u>t  +  ?  At2{at+at     (117» 


Linear  Acceleration 

«W  "  {i}t  +  J  &t  {U>t  +  5"  At'  K>t+At  tll8) 

<UW  -  {u}t +  At  tot + 1-  ^  tat  - 1-  At2  «}t+it     (119) 

Using  the  following  definitions  for  the  incremental  vectors 


Uu}t  -  {u}t+it  -  K}t  (120, 

{Au}t  =  {i}t+At  -  «}t  (121) 

{Au}t  =  {u)t+it  -  {u}t  '  (122) 


the  incremental  velocity  and  acceleration  vectors  can  be  expressed 
in  the  form 

Constant  Acceleration 

{All)  =  —  {Au>  -  {a}.  (123) 

t    .    At2     t      t 

(Au>t  =  ^|  {Au>t-  (B}t  (124) 


where 


<Ah  =  ~  {u>   +  2  {u>  (125) 

t   At    t        t 

{B>t  =  2  {ti>t  (126) 
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Linear  Acceleration 

{Au>  =  —  {Au>  -  {A>  (127) 

fc   At2     t      t 

'      {A^}t  =  At  {Au}t  "  {B}t  (128) 

where 

Mt   -  jf  «>t  ♦  3  {5}t  (129) 

{B>t  =  3  {u)t  +  j  At  {u)t  (130) 

Using  these  relations,  the  incremental  equation  of  motion,  Eq.  (115) , 
can  be  written  in  the  form 


[K]   {Au}   =  {AR}  (131) 


where 


{Ar}   =  {AR>   +  [M]  {a>   +  C  ■  [M]  {b>  (132) 

t         t  t      2  t 


[Kl   =  >[Kl   +  C   [M]  (133) 

t       t     i 


and  the  actual  incremental  displacement  vector  can  be  expressed  as 


{Au>  =  C   ({Au>  +  3{B>  )  (134) 

t     3        t         t 


Constants  C  ,  C  ,  and  C  are  given  by  the  relations 


Constant  Acceleration 


C  =  2  a  At  +  4  (135) 

1  At2  +  2  3  At 

c  =  a  -  c  B  (136) 

2  1 
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c,  -  at4V  <i37> 


Linear  Acceleration 

C  -  3  a  At  +  6  (138) 

1  At2  +  3  3  At 

c  =  a  -  c  3  (139) 

2  1 

At 

(140) 


3   At  +  33   ■ 

After  computing  the  incremental  displacement  vector  using  Eq.  (134) , 
the  corresponding  incremental  acceleration  and  velocity  vectors  are 
determined  using  Eqs.  (123)  and  (124),  or  Eqs.  (127)  and  (128),  respec- 
tively.  The  displacement,  velocity  and  acceleration  vectors  at  time 
t+At  are  then  evaluated  using  Eqs.  (120),  (121)  and  (122),  respectively. 
The  tangent  stiffness,  strain,  and  stress  for  each  element  can  now  be 
calculated  for  time  t+At. 

An  alternative  solution  of  the  incremental  equilibrium  equation, 

Eq.  (115)  can  be  obtained  by  separating  the  tangent  stiffness  matrix 

E         P 
IK]   into  its  elastic  and  plastic  parts,  [K  ]  and  [K  ]  ,  as  described 

in  Chapter  IV.   That  term  associated  with  plastic  deformation  is  then 

transferred  to  the  right  hand  side  of  the  equation  of  motion  and  is 

treated  as  an  equivalent  load  vector  [21] . 

In  the  above  described  step-by-step  integration  procedures,  the 

initial  displacements  and  velocities  at  time  t  =  0  are  assumed  equal 

to  zero.   The  very  first  incremental  acceleration  vector  is  then 

computed  directly  from  Eq.  (115),  i.e., 

•     [M]    lAu>+    [C]„{0}+ [K]      (0}        =      {AR> 
u  o  u  o  0  0 
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or 

{Au}   =   [M]"1  {AR}  =   [M]"1  {R}   =   -{u  } 
o  o  o      go 

The  initial  forces  existing  in  the  elements  cannot  be  assumed  zero  but 
must  be  taken  equal  to  the  static  gravitational  forces  since  tangent 
stiffness  is  dependent  upon  total  force  (gravity  plus  seismic) . 

While  step-by-step  integration  procedures  similar  to  those  described 
above  must  be  used  for  nonlinear  analyses/  they  may  or  may  not  be  used 
for  linear  analyses  since  the  mode  superposition  method  is  an  alternate 
method  which  can  be  used  for  linear  analyses  [116] .   In  the  present  in- 
vestigation, it  was  found  computationally  convenient  to  use  the  step- 
by-step  method  for  both  linear  and  nonlinear  analyses. 

F.   TIME  INTERVAL  At 

The  step-by-step  integration  method  is  accurate  only  if  the  time 
interval  At  is  small  compared  with  the  shortest  period  T  of  the  soil- 
structural  system  and  is  also  small  compared  with  the  predominant  pe- 
riods in  the  excitation.  Assuming  the  latter  condition  is  satisfied, 
the  ratio  At/T  must  be  selected  less  than  a  certain  critical  value  to 
insure  a  convergent  and  stable  solution  in  the  case  of  the  linear  acce- 
leration method;  however,  it  can  be  shown  that  the  constant  acceleration 
method  is  always  stable  for  a  linear  system  [76] .   In  the  present  study, 
the  presence  of  the  nonlinear  friction  elements  tend  to  encourage  an  un- 
stable response,  if  the  At/T  ratio  is  taken  too  large.   Therefore,  ex- 
treme care  must  be  taken  in  selecting  the  numerical  value  of  this  ratio. 
The  effects  of  this  particular  parameter  on  dynamic  response  are  discussed 
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subsequently  in  Chapter  VI. 

G.   EARTHQUAKE  INPUT 

In  the  present  investigation,  the  horizontal  ground  motion  was 
prescribed  in  accordance  with  the  acceleration  time-history  shown  in 
Fig.  30.  This  artificial  accelerogram  was  generated  by  A.  K.  Chopra  to 
simulate  the  ground  motions  produced  by  the  San  Fernando  earthquake  at 
the  site  of  the  Olive  View  Hospital  located  about  6  miles  southwest  of 
the  epicenter  118] .   It  has  a  peak  acceleration  of  0.5  g  and  a  uniform 
phase  of  high  intensity  shaking  for  8'  seconds. 

The  vertical  ground  motions  were  assumed  zero  for  the  present  study , 
but  the  computer  program  has  the  option  to  permit  input  of  vertical 
ground  motions. 
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VI   PARAMETER  STUDIES 

The  previously  defined  mathematical  modelling  and  dynamic  analysis 
procedures  have  been  applied  to  a  straight  version  of  an  existing  slightly- 
curved  skewed  bridge,  namely,  the  North  Connector  Undercrossing  located 
approximately  800  ft.  northerly  of  Route  5  -  San  Fernando  Road  Interchange 
in  the  city  and  county  of  Los  Angeles.   Plan  and  elvation  views  of  this 
bridge  are  shown  in  Fig.  31  along  with  cross-sectional  views  of  the  decks 
and  the  centrally  located  supporting  columns.   Figure  32  shows  a  sectional 
view  of  the  left  abutment  with  a  portion  of  the  backfill  of 'extent  L  and 
height  H. 

It  is  the  purpose  of  this  chapter  to  present  the  results  of  dynamic 
analyses  for  the  North  Connector  Undercrossing  when  subjected  to  one  com- 
ponent of  earthquake  excitation  in  its  longitudinal  direction.   The  ground 
motions  used  in  this  study  were  the  Olive  View  Hospital  accelerations 
shown  in  Fig.  30  during  the  time  interval  2-4  seconds.   This  relatively 
short  duration  was  chosen  to  minimize  computer  time  and  yet  to  provide 
sufficient  time  for  a  representative  dynamic  response  to  occur.   The 
peak  acceleration  in  this  excitation  is  approximately  0.3  g. 

In  order  to  establish  an  appropriate  mathematical  model  of  the  en- 
tire bridge-soil  system,  parameter  studies  were  first  carried  out  using 
a  rigid  wall  with  uniform  elastic  backfill.   The  complete  bridge-soil 
system  was  then  analyzed  with  certain  parameter  variations.   The  results 
of  these  studies  along  with  the  results  of  an  integration-time-interval 
sensitivity  analysis  are  presented  in  the  subsequent  sections  of  this 
chapter . 
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A.   RIGID  WALL-BACKFILL  SYSTEM 

1.   Lateral  Extent  of  Backfills  -  To  study  the  longitudinal  dimen- 
sion or  extent  of  backfill  required  in  the  mathematical  modelling  of 
bridge-soil  systems,  analyses  were  conducted  using  a  single  rigid  wall 
and  a  uniform  linear -elastic  backfill  having  the  properties  y   -  110  pcf , 
V  =  0.35,  and  G  =  10  ksi.   Five  cases  with  L/H  ratios  equal  to  4.8,  6.0, 
10.0,  14.0,  and  22.0,  as  shown  in  Fig.  33,  were  used  in  this  investigation, 
The  finite  element  idealization  consists  of  5  elements  in  the  vertical 
direction  and  a  variable  number  in  the  horizontal  direction.  The  ele- 
ment length  to  height  ratio  R  =  a/b  was  maintained  at  a  constant  value 
equal  to  2  throughout  the  extent  of  the  backfill  in  Cases  1  and  2  but 
was  maintained  only  over  a  distance  2H  from  the  wall  in  Cases  3,  4,  and 
5.   This  ratio  was  set  equal  to  4  for  all  elements  beyond  the  distance 
2H  in  the  latter  cases.  No  friction  elements  were  placed  between  the 
rigid  wall  and  the  backfill  in  these  particular  studies  and  the  left 
vertical  boundary  of  the  backfill  was  assumed  free  in  each  case.   The 
fundamental  period  of  the  soil  system  under  these  conditions  is  about 
0.08  seconds  in  each  case  with  slight  increases  occurring  with  increas- 
ing values  of  L/H  as  shown  in  Table  1. 

Assuming  rigid  body  earthquake  excitations  to  occur  along  the  en- 
tire base  of  the  backfill  and  at  the  rigid-wall  vertical-boundary  and 
assuming  5  percent  of  critical  damping  in  the  first  two  modes  of  vibra- 
tion, time  histories  of  response  were  obtained  for  all  5  cases.  Response 
quantities  of  greatest  interest  were  the  maximum  or  peak  values  of  (1) 
total  lateral  dynamic  force  exerted  on  the  rigid  wall,  and  (2)  horizontal 
displacements  and  accelerations  at  various  location  in  the  backfill. 
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Figure  34  shows  the  maximum  total  lateral  dynamic  force  F,  exerted 

d 

on  the  wall  for  each  of  the  5  cases  studied.  This  force  is  reasonably 

constant  at  a  value  of  approximately  14.5  kips  over  the  range  studied; 

i.e.  4.8  <  L/H  <  22,  and  is  about  3.9  times  greater  than  the  average 

static  value  which  is  approximately  3.7  kips.  The  mean  location  of  the 

resultant  force  F,  was  found  to  be  at  a  distance  0.52  H  above  the  base 
a 

which  is  considerably  higher  than  the  position  of  the  resultant  static 

force  F  located  at  0.34  H. 
s 

The  maximum  horizontal  displacements  relative  to  the  moving  base 
at  6  different  locations  in  the  backfill  are  shown  in  Fig.  35  for  all 
5  cases.  While  these  displacements  are  reasonably  constant  over  the 
range  of  L/H  studied,  large  differences  are  noted  from  one  point  to 
another.   The  displacements  for  points  1  and  2  are  very  small  due  to  the 
fact  they  are  located  only  one  element  away  from  the  rigid  wall.  Points 

3  and  4  which  are  5  elements  away  from  the  wall  experienced  much  larger 
displacements  than  points  1  and  2  and  points  5  and  6  at  the  left  bounda- 
ry experienced  even  larger  displacements.   These  relative  displacements 
reflect  the  manner  in  which  the  rigid  wall  boundary  effects  decay  with 
increasing  distance  from  the  wall.  The  displacements  for  points  3  and  5 
are  considerably  larger  than  the  corresponding  displacements  for  points 

4  and  6,  respectively;  thus,  demonstrating  the  increase  in  displacements 
with  vertical  distance  above  the  rigid  base. 

Figure  36  shows  the  maximum  total  horizontal  acceleration  at  loca- 
tions 1-6  for  all  5  cases  studied.   The  relative  variations  of  this 
response  quantity  with  L/H,  horizontal  distance  from  the  wall,  and  ver- 
tical distance  above  the  base  are  similar  to  those  previously  described 
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for  horizontal  displacement.   It  should  be  noted  that  the  average  accele- 
ration with  L/H  ranges  from  about  0.5  g  for  points  1  and  2  to  about  1.1  g 
for  point  5.   These  values  represent  a  rather  small  amplification  of 
acceleration  at  points  1  and  2  and  a  fairly  large  amplification  of  the 
acceleration  at  point  5  over  the  peak  excitation  acceleration  of  about 
0.3  g. 

Based  on  the  numerical  results  shown  in  Figs.  34-36  and  summarized 
in  Table  2,  it  appears  that  an  L/H  ratio  equal  to  14  is  sufficient  for 
use  in  any  bridge-soil  analysis.   This  ratio  may  even  be  reduced  to  a 
value  of  6  with  little  loss  in  accuracy  in  determining  the  maximum  values 
of  abutment  backfill  forces.  However,  this  reduction  could  introduce 
significant  changes  in  the  predicted  maximum  horizontal  accelerations. 

2.  Length  to  Height  Ratio  of  Finite  Elements  -  Length  to  height 
ratios  R  =  a/b  for  the  finite  elements  of  the  backfill  model  were  rather 
arbitrarily  assigned  values  as  shown  in  Table  1  for  the  5  cases  pre- 
viously defined.   To  investigate  the  influence  of  changing  these  ratios, 
4  new  cases  as  shown  in  Fig.  37,  each  having  an  L/H  ratio  equal  to  10, 
are  defined.   The  length  to  height  ratio  equals  2  throughout  the  model 
for  Case  1  and  over  a  distance  2H  from  the  wall  for  Cases  2,  3,  and  4. 
Beyond  2H,  the  length  to  height  ratio  equals  4,  6,  and  10  for  Cases  2, 
3,  and  4,  respectively.  Five  elements  are  again  used  in  the  vertical 
direction  of  the  model  and  fixed  and  free  boundaries  are  prescribed  at 
the  right  and  left  ends,  respectively.   Linear  elastic  soil  properties 
are  again  assumed  equal  to  the  values  previously  assigned. 

As  before,  rigid  body  earthquake  excitations  are  assumed  to  occur 
along  the  entire  base  of  the  backfill  and  5  percent  of  critical  damping 
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is  assigned  to  the  first  two  modes  of  vibration.  Time  histories  of 
response  are  obtained  for  all  4  cases,  including  (1)  total  dynamic  force 
exerted  on  the  rigid  wall,  and  (2)  horizontal  displacements  and  accelera- 
tions at  6  locations  in  the  backfill. 

Figure  38  shows  total  static  force  F  and  the  maximum  total  dynamic 

lateral  force  F,  exerted  on  the  wall  for  each  of  the  4  cases  defined 
a 

above.  This  force  is  fairly  constant  at  a  value  of  approximately  14.5 
kips  over  the  range  of  R  assigned  beyond  x  =  2H  and  is  about  3.9  times 
greater  than  the  average  static*  value.  The  mean  location  of  the  resul- 
tant force  F,  is  at  a  distance  about  0.52H  above  the  base.  The  average 
a 

position  of  the  static  resultant  is  0.34H. 

The  maximum  horizontal  displacements  relative  to  the  moving  base 
and' the  maximum  total  horizontal  accelerations  at  the  6  different  loca- 
tions are  shown  in  Figs.  39  and  40,  respectively,  for  all  4  cases.   The 
results  in  these  figures  are  quite  similar  to  the  results  shown  the 
corresponding  Figs.  35  and  36.   Thus,  it  is  apparent  that  the  changes 
introduced  in  R  for  x  >  2H  have  introduced  relatively  small  changes 
in  overall  response  and  that  Case  4,  Fig.  37,  can  be  considered  a 
reasonable  model  of  the  backfill. 

Table  3  presents  a  summary  of  maximum  response  for  the  above  des- 
cribed 4  cases. 

3.  Number  of  Finite  Elements  in  the  Vertical  Direction  -  All  pre- 
vious cases  analyzed  have  used  5  finite  elements  in  the  vertical  direc- 
tion of  the  backfill.  To  check  the  adequacy  of  this  number  the  distribu- 
tions and  resultant  magnitudes  of  the  static  and  dynamic  backfill  pres- 
sures on  the  rigid  retaining  wall  are  compared  using  5  and  10  elements 
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in  this  direction. 

Two  cases  as  defined  in  Fig.  41  are  used  for  this  comparison.   Both 
cases  use  an  L/H  ratio  equal  to  10  and  R  ratios  equal  to  "2  and  10  for 
x  <  2H  and  x  >  2H,  respectively.   Case  1  uses  5  elements  in  the  vertical 
direction  while  Case  2  uses  10  elements. 

tIt  is  quite  clear  from  the  results  in  Fig.  41  that  the  distributions 
and  maximum  resultant  magnitudes  of  the'  backfill  forces  are  very  similar 
for  Cases  1  and  2  and  that  the  positions  of  the  corresponding  resultants 
almost  coincide.  Considering  this  fact  and  the  fact  that  their  time 
histories  (see  Fig.  42)  are  very  similar,  one  may  conclude  that  5  finite 
elements  in  the  vertical  direction  of  the  backfill  is  sufficient  for 
engineering  purposes. 

4.   Soil  Stiffness  -  To  study  the  influence  of  soil  stiffness  on 
the  dynamic  response  of  the  backfill  soil  system,  the  finite  element 
model  identified  as  Case  4  in  Fig.  37  was  analyzed  for  three  soil  con- 
ditions, namely,  a  uniform  soil  modulus  equal  10  ksi,  a  uniform  soil 
modulus  equal  to  2.5  ksi,  and  a  variable  soil  modulus  in  accordance 
with  Eq.  (38)  for  K  =50.   Poisson's  ratio  V  was  assumed  equal  to  0.35 
and  the  unit  weight  was  again  assigned  the  value  110  pcf .  These  new 
soil  conditions  are  identified  as  Cases  1  -  3  in  Fig.  43. 

While  the  shapes  of  the  distributions  of  maximum  total  wall  pressures 

and  their  resultant  force  positions  are  quite  similar  for  all  three  cases, 

the  magnitudes  of  the  resultant  dynamic  wall  pressures  vary  considerably 

with  the  soil  stiffness  condition.   This  variation  ranges  from  F,= 

(100)  (13.06)  for  Case  3  to  F  =  (154.0)  (13.06)  for  Case  2.   Obviously, 

a 

for  the  particular  excitation  used,  the  less  stiff  backfill  soils  produce 
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higher  backfill  forces  on  the  rigid  wall.   If  an  excitation  having  a 
different  amplitude  distribution  with  frequency  had  heen   used,  this  ob- 
servation could  be  changed  considerably.   Therefore,  one  must  use  cau- 
tion in  interpreting  the  results  of  this  particular  study.   It  is  impor- 
tant to  recognize  however  that  soil  stiffness  can  be  an  important  para- 
meter which  should  be  studied  using  realistic  earthquake  excitations. 

5.   Frictional  Element  Stiffnesses  -  As  pointed  out  previously,  the 
shear  and  normal  stiffnesses  of  the  frictional  elements  located  between 
backfill  soil  and  abutment  walls  are  arbitrarily  assigned  finite  but 
very  high  values  rather  than  infinite  values  to  avoid  discontinuities 
in  their  force  displacement  relations.   To  study  the  influences  of  these 
stiffnesses  on  dynamic  response,  the  rigid  wall-soil  system  shown  in 
Fig.  44  was  analyzed  assuming  linear  soil  behavior.   In  these  studies 
the  shear  and  normal  stiffnesses  (K  and  K  )  for  each  frictional  element 
were  assigned  equal  values  ranging  from  1  to  10  ksi. 

The  total  static  and  total  maximum  dynamic  lateral  wall  forces 
(F  and  F  )  obtained  in  these  studies  are  plotted  in  Fig.  45  for  stiff- 
nesses K  =  K_  equal  to  1,  10  ,  10b,  and  10  ksi.  As  one  would  expect, 
these  forces  are  reasonably  constant  over  a  wide  range  of  stiffnesses. 

However,  as  the  stiffnesses  approach  zero,  F  and  F,  also  approach  zero 

s      d 

which  represent  unrealistic  values.  Theoretically  both  F  and  F,  should 
approach  constant  values  asymptotically  with  increasing  stiffnesses, 
however  the  value  of  F  for  K  =  K  =  10  ksi  is  much  larger  than  for 
the  smaller  values  of  stiffness.  This  large  increase  is  due  to  a  numeri- 
cal instability  which  developed  in  the  analysis  procedures  and  therefore 
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should  be  ignored.   It  appears  therefore  that  realistic  wall  forces  can 
be  obtained  by  selecting  stiffnesses  in  the  range  10  <  K  =  K  <  10 
ksi. 

Although  the  asymptotic  static  value  of  wall  force  in  Fig.  45  is 
consistent  with  values  previously  presented  for  cases  having  no  friction 
elements,  the  maximum  dynamic  force  of  about  40.0  kips  is  considerably 
larger  than  the  average  value  (14.5)  previously  presented.  This  increase 
in  dynamic  wall  force  is  due  to  the  separations  and  associated  impacts 
which  occur  between  the  backfill  soil  and  the  upper  part  of  the  rigid 
wall.   Thus  it  appears  that  for  high  intensity  excitations,  the  friction 
element  is  essential  to  realistic  modelling. 

B.   INTEGRATION-TIME-INTERVAL  SENSITIVITY  ANALYSIS 

Throughout  the  rigid  wall-backfill  parameter  studies  previously 
described,  the  numerical  integration  time  step  was  assigned  a  value 
equal  to  0.01  seconds  and  the  constant  acceleration  method  was  used  which 
enables  response  to  be  stable,  but  not  necessarily  convergent,  for  all 
modes  of  vibration.   To  study  the  adequacy  of  using  0.01  seconds  for  At, 
the  rigid  wall-soil  system  defined  by  Case  4,  Fig.  37,  was  re-analyzed 
using  At  =  0.001  seconds.   The  total  number  of  degrees  of  freedom  for 
this  system  is  90  with  the  fundamental  period  being  0.084  seconds  and 
the  highest  period  estimated  at  0.0065  seconds.   The  convergent  limit 
of  the  ratio  of  time  step  duration  to  period,  i.e.  At/T,  is  0.39  for  the 
constant  acceleration  method.   Thus,  for  At  =  0.01  seconds,  the  conver- 
gent period  T  is  0.026  seconds.   Since  this  period  corresponds  to  the 
period  of  the  22nd  mode  of  vibration,  only  the  lowest  22  modes  are 
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convergent  in  the  constant  acceleration  method  of  analysis  when  At  =  0,01 
seconds.   If  on  the  other  hand,  At  =  0.001  seconds,  all  modes  are  con- 
vergent . 

To  check  the  accuracy  of  response  obtained  for  the  above  case  using 
At  =  0.01  seconds,  the  time  histories  of  lateral  dynamic  wall  force  and 
horizontal  acceleration  at  point  No.  1  (see  Fig.  40)  are  obtained  for 
At  =  0.001  and  are  plotted  in  Fig.  46  where  they  can  be  compared  with 
corresponding  results  for  At  =  0.01  seconds.  The  vertical  acceleration 
time  histories  for  point  No.  1  are  plotted  in  Fig.  47  for  At  =  0.01 
and  0.001  seconds.  A  summary  of  the  maximum  values  of  response  for 
this  case  is  presented  in  Table  4. 

Obviously,  the  results  of  Fig.  46  indicate  that  a  time  step  inter- 
val of  0.01  seconds  is  quite  adequate  in  predicting  total  lateral  wall 
force  and  horizontal  acceleration  time  histories.  However,  the  results 
of  Fig.  47  indicate  the  very  low  level  vertical  acceleration  time  his- 
tories caused  primarily  by  very  high  frequency  modal  responses  cannot 
be  predicted  accurately  by  At  =  0.01  seconds.   Since  this  high  frequency 
response  is  relatively  unimportant  from  an  engineering  point  of  view, 
it  is  concluded  that  At  =  0.01  seconds  is  adequate  for  the  previously 
described  rigid  wall-soil  parameter  studies  and  also  for  the  bridge-soil 
system  studies  to  be  described  subsequently. 

C.   BRIDGE-SOIL  SYSTEM 

To  study  the  dynamic  response  of  the  combined  bridge-soil  system, 
3  mathematical  models  of  the  North  Connector  Undercrossing  as  shown  in 
Fig.  48  were  defined.   Model  A  has  fixed  boundary  conditions  at  depth  H 

-85- 


of  the  backfilis,  at  the  base  of  abutments,  and  at  the  base  of  all  col- 
umns. Model  B  has  fixed  boundary  conditions  at  depths  2.2H  and  2.5H  of 
the  backfills,  leveling  with  bases  of  pier  columns* which  allows  the  base 
of  abutments  to  translate  and  rotate  with  the  soil  system,  and  fixed 
boundary  conditions  are  also  provided  at  the  base  of  all  columns.  Model 
C  has  fixed  boundary  conditions  only  along  the  base  of  the  backfills  as 
in  Model  B.  The  bases  of  abutments  and  columns  of  this  model  are  attached 
to  equivalent  columns  representing  the  foundation  flexibility.  These 
equivalent  columns  of  course,  have  fixed  boundary  conditions  at  their 
bases. 

The  soil  elements  in  all  three  models  can  be  assumed  linear  or  non- 
linear as  desired  and  friction  elements  can  be  included  in  Models  A  and 
C,  but  not  in  Model  B.  The  backfills  extend  a  distance  6H  in  all  models 
as  shown.  The  bridge  deck  is  linear  in  each  model;  however,  either 
linear  or  nonlinear  columns  can  be  used.   Backfill  and  foundation  soils 
were  assumed  to  have  the  properties  G  =  10.0  ksi,  V  =  0.35,  y   =  110 
lb/ft3,  c  =  0,  and  <f>  =  30°. 

1.   Soil  Pressures  on  Abutments  -  The  static  and  maximum  dynamic 

pressure  distributions  on  one  abutment  wall  for  two  cases  are  shown  in 

Fig.  49.  The  model  used  for  Case  1  is  Model  A  with  no  friction  elements 

and  with  all  other  elements  assumed  linear.  The  model  used  for  Case  2 

is  Model  B  in  its  complete  linear  form.  Due  to  the  characteristic 

response  of  the  bridge-soil  system,  the  static  and  dynamic  pressure 

distributions  are  quite  different  in  form  in  each  case.   The  resultant 

lateral  static  force  F  for  Case  2  is  about  17%  less  than  for  Case  1, 

s 

due  to  the  change  in  abutment  flexibility  and  the  maximum  resultant 
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dynamic  force  F  for  Case  2  exceeds  the  value  for  Case  1  by  170%.  This 
latter  difference  is  undoubtedly  due  to  a  closer  matching  of  the  lower 
mode  periods  of  vibration  for  Case  2  with  the  predominant  periods  in  the 
excitation.  The  location  of  the  resultant  static  force  is  at  about  H/3 
for  Case  1  and  at  about  H/2  for  Case  2  while  the  location  of  the  maxi- 
mum dynamic  resultant  force  is  at  about  0.53H  and  0.6H,  respectively. 

Another  check  on  the  influence  of  bridge  structure  flexibility  on 
the  resultant  abutment  soil  forces  can  be  made  by  comparing  the  results 
for  Cases  1  and  2  in  Fig.  50.   In  this  figure.  Case  1  is  identical  to 
Case  1  in  Fig.  49  but  Case  2  is  different.  Here,  Case  2  is  actually  the 
same  as  Case  1  except  that  the  abutment  wall  and  column  stiffnesses 
have  been  reduced  by  a  factor  of  10.   It  is  seen  that  Case  2  shows  a  75% 
increase  in  the  maximum  dynamic  resultant  force  over  Case  1  due  to  the 
decrease  in  bridge  structure  flexibility.  This  increase  is  consistent 
with  the  similar  increase  previously  noted  for  Case  2  in  Fig.  49.  The 
location  of  the  resultant  dynamic  lateral  force  is  again  at  about  mid- 
height.  The  increase  in  the  resultant  static  force  for  Case  2  over  the 
value  for  Case  1  is  due  to  the  increase  in  rotation  (due  to  deck  dead 
loads)  at  the  top  of  the  abutment  caused  by  the  reduced  abutment  flex- 
ibility. 

Further  results  of  analysis  are  shown  in  Fig.  51  identified  as  Cases 
1,  2,  and  3.   In  this  figure,  Case  1  is  again  the  complete  linear  ver- 
sion of  Model  A.  Cases  2  and  3  are  also  based  on  Model  A  but  Case  2 
has  introduced  one  nonlinearity,  namely  the  friction  elements,  and  Case 
3  has  employed  two  nonlinearities  -  friction  elements  and  nonlinear 
soil  elements.  The  very  large  increase  in  maximum  dynamic  force  F  for 
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Cases  2  and  3  over  Case  1  is  due  primarily  to  the  impact  wall  forces 
following  separations  between  wall  and  backfill. 

Finally  the  results  of  two  additional  analyses  are  s,hown  in  Figure 
52.   Cases  1  and  2  in  this  figure  are  based  on  Model  C  using  friction 
elements  and  linear  soil  elements;  however,  Case  2  uses  the  equivalent 
foundation  columns  while  Case  1  does  not.  The  most  significant  result 
to  note  in  Fig.  52  is  that  force  F  is  more  than  twice  as  great  for  Case 
2  over  the  value  shown  for  Case  1.  Again  this  increase  is  due  to  the 
fact  that  the  more  flexible  bridge  system  has  a  closer  matching  of  fre- 
quencies with  the  predominant  frequencies  in  the  excitation. 

2.   Total  Seismic  Force  Carried  by  Columns  and  Abutments  -  To  inves- 
tigate the  maximum  total  base  shear  carried  by  columns  and  abutments, 
five  cases  were  analyzed  as  indicated  in  Table  5.  Case  1  represents  the 
bridge  alone  with  no  soil-structure  interaction,  i.e.  Model  A,  Fig.  48, 
but  with  no  backfill;  Case  2  is  Model  A  with  linear  backfill  and  no  fric- 
tion elements;  Case  3  is  Model  .B*  Fig.  48,  with  linear  backfill;  Case  4 
is  Model  A  with  linear  backfill  and  friction  elements;  and  Case  5  is 
Model  A  with  non-linear  backfill  and  friction  elements.  Clearly,  the 
presence  of  backfill  contributes  significantly  to  the  maximum  total  base 
shear,  also  it  appears  from  Case  3  that  the  total  base  shear  increases 
with  overall  flexibility  which  is  again  evidence  of  a  better  matching  of 
the  lower  natural  frequencies  with  the  predominant  frequencies  in  the 
excitation.  The  relatively  large  displacement  shown  for  Case  3  is  due  to 
the  large  flexibility  of  the  system  for  this  case  in  comparison  with  the 
other  cases. 
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Do   COMPARISON  OF  RESULTANT  ABUTMENT  BACKFILL  FORCE  OBTAINED  BY  ANALYSIS 
AND  THE  MONONOBE-OKABE  METHOD 


One  commonly  used  formula  in  calculating  the  resultant  dynamic  lateral 
force  on  the  abutment  wall  is  the  Mononobe-Okabe  formula  [53,  55,  97] . 
This  formula  has  the  following  form 

Pp  =  {1-Kv)  *  r  '  x  *  V  (141) 

where 

P  =  Passive  earthpressure  at  depth  x 

K  =  vertical  seismic  coefficient 
v 

Y  =  Unit  weight  of  soil 

x  =  arbitrary  depth 

K_  =  Passive  earthpressure  coefficient  during  earthquake 

The  coefficient  K   is  given  by 

.     cos2    (<J>-e0+8) 


Kg  r  /sin  $    •    sin(<f>+a-U0)  "| 2 

cos(6-a) 


[/sin  <p   •    sin 
1  -  / ,Q  Q   . — 
cos  (6-0  0)    • 


(142) 


where 


<J>  =  Angle  of  friction  of  soil 
—  l   ^h 

eo  =  tan   I3T 

v 

1L  =  Horizontal  seismic  coefficient 

0  =  Angle  between  the  backline  of  the  wall  and  the  vertical 
line 
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a  =  Angle  between  ground  surface  and  the  horizontal  line 

'o  -  tan   IT^o 


Using  <p  =   30°;  0  =  0°;  9  =  tan  *  ,°'1n  =   6°;  a  =  0°;  y  =   110  lb/ft3 


K   =2.86 

The  earthpressure  at  bottom  of  wall  is 

p  =  (1-0.0)  *  110  *  13.5  *  2.86  =  4.25  K/ft2 
P 

and  the  total  lateral  force  on  the  wall  is 
FJ  =  i  *  13.5  *  4.25  =  28.6  K/ft  ' 

Using  the  non-linear  form  of  Model  A,  Fig.  48,  or  Case  3  of  Fig.  51 , 
analysis  gives 

P.  =  35.0  K/ft 

a 

■This  analytical  result  is  a  higher  value  than  that  formula  given  by  the 
Monobe-Okabe .  The  position  of  the  resultant  force  is  assumed  a  dis- 
tance H/3  above  the  base  when  using  the  formula;  however,  the  analysis 
shows  it  to  be  0.44H  above  the  base.   This  higher  position  given  by  an 
analysis  is  consistent  with  other  investigations  [97,  117]. 
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VII     GENERAL  CONCLUSION 

Based  on  the  results  of  this  investigation,  it  is  concluded  that 
soil- structure  interaction  effects  must  be  considered  when  analyzing  the 
dynamic  response  of  short,  stiff,  single  or  multiple  span  bridges.  The 
mathematical  modelling  and  computer  programs  presented  herein  provide  an 
effective  means  of  conducting  such  analyses. 

Since  the  numerical  results  obtained  in  this  investigation  are  very 
limited,  caution  should  be  exercised  when  interpreting  them  in  a  quanti- 
tative sense.  Further  analyses  are  recommended  to  complete  the  parameter 
studies. 
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